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ABSTRACT 


Title  of  Dissertation:  APPROXIMATE  MATRIX  DIAGONALIZATION 

FOR  USE  IN  DISTRIBUTED  CONTROL 
NETWORKS 

George  A.  Kantor,  Doctor  of  Philosophy,  1999 

Dissertation  directed  by:  Professor  P.S.  Krishnaprasad 

Department  of  Electrical  and  Computer  Engineering 


Distributed  control  networks  are  rapidly  emerging  as  a  viable  and  important 
alternative  to  centralized  control.  In  a  typical  distributed  control  network,  a  num¬ 
ber  of  spatially  distributed  nodes  composed  of  “smart”  sensors  and  actuators  are 
used  to  take  measurements  and  apply  control  inputs  to  some  physical  plant.  The 
nodes  have  local  processing  power  and  the  ability  to  communicate  with  the  other 
nodes  via  a  network.  The  challenge  is  to  compute  and  implement  a  feedback  law 
for  the  resulting  MIMO  system  in  a  distributed  manner  on  the  network. 

Our  approach  to  this  problem  is  based  on  plant  diagonalization.  To  do  this,  we 
search  for  basis  transformations  for  the  vector  of  outputs  coming  from  the  sensors 


and  the  vector  of  inputs  applied  to  the  actuators  so  that,  in  the  new  bases,  the 
MIMO  system  becomes  a  collection  of  decoupled  SISO  systems.  This  formulation 
provides  a  number  of  advantages  for  the  synthesis  and  implementation  of  a  feed¬ 
back  control  law,  particularly  for  systems  where  the  number  of  inputs  and  outputs 
is  large.  Of  course,  in  order  for  this  idea  to  be  feasible,  the  required  basis  transfor¬ 
mations  must  have  properties  which  allow  them  to  be  implemented  on  a  distributed 
control  network.  Namely,  they  must  be  computed  in  a  distributed  manner  which 
respects  the  spatial  distribution  of  the  data  (to  reduce  communication  overhead) 
and  takes  advantage  of  the  massive  parallel  processing  capability  of  the  network 
(to  reduce  computation  time). 

In  this  thesis,  we  present  some  tools  which  can  be  used  to  find  suitable  trans¬ 
forms  which  achieve  “approximate”  plant  diagonalization.  We  begin  by  showing 
how  to  search  the  large  collection  of  orthogonal  transforms  which  are  contained 
in  the  wavelet  packet  to  find  the  one  which  most  nearly,  or  approximately,  diag¬ 
onalizes  a  given  real  valued  matrix.  Wavelet  packet  transforms  admit  a  natural 
distributed  implementation,  making  them  suitable  for  use  on  a  control  network. 
We  then  introduce  a  class  of  linear  operators  called  recursive  orthogonal  transforms 
(ROTs)  which  we  have  developed  specifically  for  the  purpose  of  signal  processing 
on  distributed  control  networks.  We  show  how  to  use  ROTs  to  approximately 
diagonalize  fixed  real  and  complex  matrices  as  well  as  transfer  function  matrices 
which  exhibit  a  spatial  invariance  property.  Numerical  examples  of  all  proposed 
diagonalization  methods  are  presented  and  discussed. 
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Chapter  1 


Introduction 


Economics  often  plays  a  deciding  role  in  the  design  and  implementation  of  a  large 
control  system.  In  the  past,  communication  was  cheap  compared  to  computational 
power.  As  a  result,  centralized  control  became  the  norm,  with  one  expensive 
central  computer  reading  the  data  directly  from  all  of  the  sensors  and  sending 
control  inputs  directly  to  all  of  the  actuators. 

Today,  however,  things  are  changing.  Mass  produced  microprocessors  now 
provide  continually  increasing  amounts  of  computational  power  at  continually  de¬ 
creasing  prices.  In  automotive  applications,  for  example,  the  cost  of  installing  and 
maintaining  the  wiring  required  for  a  large  centralized  control  system  has  exceeded 
the  cost  of  the  computational  equipment  [38].  Additionally,  the  centralized  con¬ 
troller  is  difficult  to  reconfigure:  adding  a  new  sensor  or  actuator  to  the  plant 
requires  the  installation  of  a  new  set  of  wires  connecting  the  new  element  to  the 
computer.  Hence,  it  makes  sense  to  investigate  alternatives  to  centralized  control. 

Distributed  control  networks  have  emerged  as  a  viable  alternative  to  centralized 
control.  In  such  a  network,  “intelligent”  sensors  and  actuators  make  control  deci¬ 
sions  based  on  a  combination  of  local  information  and  digital  commands  from  the 
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network.  Each  node  is  also  capable  of  sending  commands  and  important  sensor  in¬ 
formation  to  the  other  nodes  on  the  network.  For  example,  an  “intelligent”  motor 
might  receive  a  position  command  over  the  network,  then  use  an  on-board  trajec¬ 
tory  planner,  sensor,  and  PID  controller  to  get  to  that  position.  Finally,  when  the 
motor  has  reached  its  desired  position,  it  would  broadcast  its  new  position  to  the 
other  sensors  and  actuators  on  the  network. 

Control  networks  have  many  advantages  over  centralized  control.  Less  wiring  is 
required  to  connect  the  elements  of  the  system  together,  making  the  network  easier 
and  cheaper  to  install.  The  simplicity  of  the  wiring  scheme  makes  maintenance 
much  easier;  with  less  wiring,  there  is  less  to  go  wrong.  The  on-board  “intelligence” 
of  the  sensors  and  actuators  allows  them  to  perform  self  diagnostics,  a  feature 
which  eases  troubleshooting.  The  system  is  easily  reconhgurable.  New  sensors  and 
actuators  can  easily  be  connected  to  the  existing  network.  Obsolete  nodes  can  be 
replaced  in  the  same  way. 

Much  work  has  been  done  regarding  control  networks  over  the  past  few  years. 
Brockett  has  investigated  the  stabilization  of  a  network  of  intelligent  motors  [9]. 
Wong  and  Brockett  have  studied  the  problems  of  state  estimation  and  feedback 
control  for  control  networks  with  limited  communication  bandwidth  [47]  [46]. 
Hristu  [24]  has  addressed  the  problem  of  Ending  stabilizing  feedback  laws  for  lin¬ 
ear  systems  with  limited  communication.  Wang  and  Mau  [41]  and  Ooi,  Verbout, 
Ludwig,  and  Wornell  [32]  present  some  results  on  the  characteristics  of  feedback 
systems  where  the  feedback  data  is  subject  to  a  communication  delay.  Many  au¬ 
thors  have  investigated  the  problem  of  modeling  and  controlling  hybrid  systems 
[19]  [5]  [6]  [7]  [8].  Most  of  this  work  is  motivated  by  the  hybrid  dynamics  that  re¬ 
sult  from  the  combination  of  continuous  physical  interaction  and  discrete  network 
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interaction  that  takes  place  in  a  control  network. 

There  is  a  broad  range  of  applications  of  control  networks  that  are  being  put 
into  use  today.  Automated  homes  and  offices  use  control  networks  to  automatically 
turn  lights  on  and  off,  control  temperature,  manage  security  systems,  and  automate 
a  number  of  other  operations  [34],  An  control  network  can  provide  a  manufacturing 
plant  with  an  assembly  line  that  is  easily  reconfigured  to  make  different  products. 
Automobile  manufacturers  employ  control  networks  to  coordinate  anti-lock  brak¬ 
ing,  engine  management,  traction  control,  powered  seats  and  windows,  pollution 
control,  and  a  variety  of  other  systems  which  reside  on  a  modern  car  [38].  Oil 
refineries  and  chemical  plants  use  control  networks  with  intelligent  instrumenta¬ 
tion  capable  of  remote  calibration  and  automated  troubleshooting,  minimizing  the 
control  engineer’s  visits  to  inhospitable  or  dangerous  areas  of  the  plant  [25]. 

Here,  we  concentrate  on  applications  of  control  networks  where  the  number 
of  inputs  and  outputs  is  large  and  the  interaction  between  them  is  linear.  Our 
interest  in  this  problem  is  motivated  by  recent  developments  in  the  fields  of  micro- 
electro-mechanical  systems  (MEMS)  and  smart  structures.  In  the  held  of  MEMS, 
significant  advances  have  been  made  in  the  construction  of  sensors  and  actuators 
on  a  very  small  scale.  For  example,  Bifano  et.al.  [3]  have  fabricated  a  MEMS  mirror 
array  for  optical  image  processing.  The  array  fits  400  actuators  onto  an  area  of 
less  than  1  square  centimeter.  Since  most  MEMS  devices  are  currently  fabricated 
in  silicon,  the  idea  of  designing  MEMS  sensors  and  actuators  with  built-in  local 
microprocessors  seems  natural.  In  the  held  of  smart  structures,  a  large  number  of 
sensors  and  actuators  are  incorporated  into  the  natural  structure  of  a  mechanical 
system.  One  such  application  has  large  number  of  piezoelectric  actuators  and 
strain  sensors  embedded  in  a  composite  panel.  The  panel  could  then  be  a  hight 
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surface  on  an  aircraft  and  the  sensor  actuator  system  would  be  used  to  control 

its  shape,  limit  vibrations,  and  monitor  its  health.  With  such  a  large  number  of 

sensors  and  actuators  involved,  the  implementation  of  a  controller  on  a  distributed 

control  network  seems  more  feasible  than  the  prospect  of  a  centralized  controller. 

A  similar  panel  could  be  placed  in  the  interior  of  an  aircraft  or  automobile  and 

the  control  network  would  be  used  to  cancel  noise  and  vibrations  coming  from 

outside.  Similar  approaches  have  been  proposed  to  damp  vibrations  and  control 

the  flight  surface  of  a  helicopter  blade.  Another  application  uses  “smart”  actuators 

TM 

composed  of  an  underwater  acoustic  actuator  and  sensor.  The  SmartPanel 
[17]  and  Composite  Smart  Material  (CSM)  Tile  [45]  are  examples  of  two  such 
devices  that  are  currently  under  development.  Thousands  of  these  “smart  tiles” 
will  be  mounted  on  the  outside  of  the  hull  of  a  submarine,  covering  the  entire  hull. 
This  massive  sensor /actuator  array  will  then  be  used  to  actively  reduce  acoustic 
radiation,  cancel  enemy  sonar  pulses,  and  perform  acoustic  sensing. 

Clearly,  there  is  a  need  to  systematically  address  problems  of  this  type.  In 
this  thesis,  we  present  some  results  to  aid  in  the  design  and  implementation  of 
feedback  controllers  for  systems  equipped  with  large  distributed  control  networks. 
The  ideas  here  are  designed  to  take  advantage  of  the  parallel  processing  capability 
of  the  network  while  reducing  the  need  for  global  communication. 

1.1  Plant  Diagonalization 

The  central  theme  of  this  thesis  is  what  we  call  plant  diagonalization.  Here  we 
define  the  structure  of  the  plants  we  will  address.  We  then  introduce  the  concept 
of  plant  diagonalization  and  discuss  the  advantages  that  this  idea  presents  in  the 
context  of  control  networks. 
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1.1.1  Plant  Structure 


In  this  thesis,  we  will  consider  plants  of  the  form 

V  =  Gu,  (1.1) 

where  y  is  an  n-dimensional  vector  containing  the  plant  outputs,  a  is  an  m-dinren- 
sional  vector  of  plant  inputs,  and  G  is  the  nxm  matrix  which  defines  the  interaction 
between  the  inputs  and  outputs.  We  call  G  the  “plant  matrix” . 

In  most  of  this  thesis,  we  will  consider  plant  matrices  whose  elements  are  con¬ 
stant  real  or  complex  numbers.  These  types  of  matrices  can  be  used  to  model  two 
types  of  systems:  “quasi-static”  systems  and  “fixed  frequency”  systems.  The  term 
“quasi-static”  is  used  to  describe  stable  systems  whose  transients  are  negligible 
compared  to  the  DC  response  of  the  system.  The  plant  matrix  of  a  quasi-static 
system  is  a  constant  real-valued  matrix.  A  “fixed  frequency”  system  is  linear  dy¬ 
namic  system  whose  inputs  (and  any  disturbances)  are  assumed  to  be  sinusoids 
of  a  constant  frequency.  The  plant  matrix  of  a  fixed  frequency  system  is  the 
complex-valued  matrix  that  results  from  evaluating  the  transfer  function  matrix 
of  the  dynamical  system  at  the  input  frequency.  The  (i,j)th  element  of  the  plant 
matrix  of  a  fixed  frequency  system  is  a  complex  number  which  represents  the  gain 
and  phase  shift  from  the  jth  input  to  the  Ah  output.  Models  of  this  type  are  useful 
for  systems  which  exhibit  strong  resonances  as  will  be  discussed  in  Section  2.4.2. 

In  addition  to  this  basic  structure,  we  assume  that  the  inputs  and  outputs  cor¬ 
respond  to  actuators  and  sensors  spatially  distributed  on  a  control  network.  The 
network  is  composed  of  a  number  of  nodes.  Each  node  contains  a  combination  of 
sensors  and  actuators  along  with  a  microprocessor  and  the  communication  hard¬ 
ware  required  to  communicate  with  other  nodes  on  the  network.  Each  node  has 
access  to  its  own  data;  it  can  read  the  outputs  of  the  sensors  it  contains  and  it  can 
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apply  inputs  to  the  actuators  it  contains.  It  shares  data  with  other  nodes  only  by 
communicating  via  the  network. 

1.1.2  Diagonalizing  the  Plant  Matrix 

The  goal  of  this  thesis  is  to  find  an  efficient  method  of  implementing  output  feed¬ 
back  control  so  that  the  closed  loop  system  achieves  some  desired  performance 
criteria.  Our  approach  to  solving  this  problem  is  to  diagonalize  the  plant  matrix 
G.  To  do  this,  we  search  for  invertible  basis  transformations  for  the  input  and 
output  vectors  so  that,  in  the  new  bases,  the  transformed  plant  matrix  appears 
diagonal.  Here,  we  call  a  matrix  E  diagonal  if  [E]^  =  0  whenever  i  ^  j.  Using  this 
terminology,  it  makes  sense  to  call  a  non-square  matrix  diagonal. 

Consider  the  basis  transformations  y  =  Qiy  and  u  =  Q2u.  In  these  new  bases, 
the  input/output  behavior  of  the  plant  is  described  by  the  equation 

y  =  Q1GQz1u.  (1.2) 

Hence,  the  problem  of  finding  suitable  basis  transformations  is  equivalent  to  finding 
Q i  and  Q2  such  that  the  matrix 

E  =  QXGQ^  (1.3) 

is  diagonal. 

Suppose  it  is  possible  to  find  suitable  Q\  and  Q2.  The  system  in  the  new 
coordinates  can  be  written 

y  =  Zu.  (1.4) 

Since  E  is  diagonal,  the  interaction  between  ut  and  y3  is  zero  for  i  ^  j.  The  plant 
has  become  a  system  of  decoupled  scalar  subsystems, 

Vi  =  (Jim,  (1.5) 
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for  i  —  1, 2, . . .  ,min(n,  m). 

This  representation  has  two  major  advantages,  particularly  when  n  and  m 
are  large.  First,  the  problem  of  MIMO  controller  synthesis  is  transformed  to  the 
much  easier  problem  of  synthesizing  a  controllers  for  a  number  of  decoupled  SISO 
systems.  The  other  advantage  is  closely  related.  Since  the  systems  are  decoupled, 
the  feedback  laws  are  also  decoupled.  Each  Tq  can  be  calculated  based  only  on  the 
value  of  the  corresponding  f/j.  As  a  result,  once  y  is  computed,  u  can  be  computed 
quickly,  in  parallel,  and  with  no  additional  communication. 

To  illustrate  this,  we  compare  the  task  of  computing  the  linear  feedback  law 
u  =  Ky  in  the  original  basis  with  the  task  of  computing  u  =  Ky  in  the  transformed 
basis.  In  the  original  basis,  u  is  computed  by  multiplying  the  n  -vector  y  by 
the  m  x  n  matrix  K.  This  calculation  requires  0(nm )  (“on  the  order  of  nm”) 
operations.  Additionally,  global  communication  is  required;  in  order  to  compute  iq, 
the  value  of  every  element  of  y  must  be  known.  In  the  transformed  basis,  however, 
the  matrix  K  is  diagonal,  so  u  can  be  computed  in  (9(min(m,  n ))  operations.  Due 
to  the  decentralized  nature  of  K ,  no  additional  communication  is  required  to  make 
the  calculations  and  the  calculations  can  be  easily  performed  in  parallel  on  the 
control  network. 

The  computational  and  communication  requirements  of  the  transforms  Q\  and 
Q 2  are  very  important.  In  addition  to  the  cost  involved  in  computing  u  from  y 
we  must  also  perform  the  coordinate  transformations  y  =  Qiy  and  u  =  Q^u.  If 
these  transforms  are  computationally  intensive  or  require  global  communication, 
then  the  advantages  of  the  transformed  basis  are  lost.  In  most  cases,  this  means 
that  we  will  have  to  settle  for  transforms  which  “approximately”  diagonalize  the 
plant  matrix. 
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In  this  thesis,  we  will  implement  the  approximations  of  Q i  and  Q 2  by  using 
a  series  of  alternating  communication  and  computation  steps.  In  the  communi¬ 
cation  steps,  a  limited  amount  of  data  is  shared  between  selected  nodes.  In  the 
computation  steps,  the  nodes  individually  perform  a  simple  manipulation  of  the 
data  they  have  access  to.  By  their  design,  these  approximations  take  advantage 
of  the  distributed  computing  power  available  on  the  network  while  respecting  the 
communication  constraints. 

Distributed  Controller  Implementation 

Now  we  are  ready  to  describe  how  plant  diagonalization  can  be  applied  to  the 
implementation  of  feedback  on  a  distributed  control  network.  In  the  proposed  im¬ 
plementation,  the  controller  performs  three  basic  tasks.  First,  it  computes  y  =  Q 1  y 
using  a  series  of  alternating  communication  and  computation  steps  as  described 
above.  The  elements  of  the  output  vector  y  are  the  measurements  of  the  sen¬ 
sors  that  are  distributed  spatially  over  the  network.  Because  of  of  the  distributed 
nature  of  our  implementation  of  the  transform  Q 1;  the  elements  of  the  resulting 
transformed  vector  y  are  also  distributed  over  the  nodes  of  the  network.  The  sec¬ 
ond  task  of  the  controller  is  to  compute  the  transformed  output  vector  u  based  011 
y.  The  plant  matrix  is  diagonal  (or  “approximately  diagonal” )  in  the  transformed 
coordinates,  so  each  Ui  can  be  chosen  based  soley  on  the  corresponding  y,.  This 
means  that  each  Ui  can  be  computed  on  the  node  on  which  yr  resides  without 
any  additional  communication.  The  third  and  final  task  of  the  controller  imple¬ 
mentation  is  to  transfrom  the  vector  u  into  the  actual  input  vector  u  using  the 
transformation  u  =  Q^u.  Like  Q 1,  Q^1  is  implemented  in  a  distributed  manner 
using  a  series  of  communication  and  computation  steps.  After  the  third  task  is 


completed,  each  element  Ui  of  the  input  vector  will  reside  on  the  node  containing 
the  actuator  to  which  iq  will  be  applied. 

To  summarize,  the  main  goal  of  this  thesis  is  to  find  linear  coordinate  trans¬ 
formations  Q  i  and  Q->  which  meet  the  following  criteria: 

1.  The  matrix  Q1GQ2  1  is  diagonal  or  “approximately”  diagonal 

2.  The  transforms  Qi  and  Q2  can  be  implemented  as  a  series  of  alternating  com¬ 
putation  and  communication  steps.  The  computations  necessary  to  perform 
the  transform  can  be  easily  distributed  over  the  processors  on  the  control 
network,  taking  advantage  of  the  network’s  parallel  processing  capability. 
The  communication  steps  can  be  chosen  to  respect  the  spatial  distribution 
of  the  data,  reducing  communication  overhead. 

1.2  Preview  of  the  Thesis 

In  the  following  pages,  we  present  some  methods  of  finding  suitable  transforms  to 
diagonalize  or  approximately  diagonalize  various  types  of  plants.  In  Chapter  2,  we 
show  how  to  approximately  diagonalize  a  constant  real  or  complex  valued  matrix 
using  transforms  from  the  wavelet  packet.  The  idea  is  inspired  by  the  work  of 
Chou,  Guthart,  and  Flamm  [14],  who  proposed  the  use  of  the  wavelet  transform 
in  the  spatial  domain  to  simplify  the  required  feedback  for  vibrating  systems.  The 
work  presented  here  uses  an  algorithm  due  to  Wickerhauser  [43]  [44]  to  find  wavelet 
packet  transforms  which  approximate  the  SVD  factors  of  the  plant  matrix.  We 
show  that  the  resulting  transforms  have  a  natural  implementation  on  a  hierarchical 
control  network. 

In  Chapter  3,  we  introduce  a  class  of  transforms  which  we  have  designed  specif- 
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ically  to  be  implemented  on  a  control  network.  We  begin  by  developing  signal  pro¬ 
cessing  schemes  to  be  implemented  on  networks  whose  nodes  are  sensor/actuator 
pairs.  We  then  generalize  these  ideas  to  create  a  class  of  transforms  suitable  for 
general  control  networks.  We  call  these  transforms  recursive  orthogonal  trans¬ 
forms,  or  ROTs.  Recursive  orthogonal  transforms  are  the  subject  of  the  rest  of  the 
thesis. 

In  Chapter  4  we  consider  the  problem  of  finding  an  ROT  which  approximately 
diagonalizes  a  fixed  real  or  complex  valued  plant  matrix.  We  first  consider  fixed 
plant  matrices  which  are  real  and  symmetric.  The  approach  we  take  is  very  sim¬ 
ilar  to  Brockett’s  work  on  diagonalizing  symmetric  matrices  via  gradient  flows  on 
orthogonal  matrices  [10].  We  then  consider  the  problem  of  diagonalizing  general 
complex  matrices.  Here  we  extend  the  work  of  Helmke  and  Moore  [21]  to  find  ROTs 
which  most  closely  approximate  the  complex  SVD  factors  of  the  plant  matrix. 

In  Chapter  5,  plant  matrices  which  exhibit  a  spatial  invariance  property  are 
considered.  This  spatial  invariance  property  is  analogous  to  the  time  invariance 
of  LTI  systems.  Spatial  invariance  can  be  exploited  to  aid  in  the  design  and 
implementation  of  feedback  controllers.  Our  work  here  is  motivated  by  a  number  of 
authors,  including  El-Sayed  and  Krishnaprasad  [16]  and,  more  recently,  Bamieh  [1], 
We  show  that  all  spatially  invariant  systems  are  exactly  diagonalized  by  the  discrete 
Fourier  transform.  Armed  with  this  knowledge,  we  seek  to  End  an  ROT  capable 
of  diagonalizing  any  spatially  invariant  plant  matrix. 

Chapter  6  contains  a  quantitative  comparison  of  the  approximate  matrix  di- 
agonalization  methods  presented  in  this  thesis.  We  also  make  some  concluding 
remarks  and  give  some  suggestions  for  future  work. 
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Chapter  2 


Approximate  SVD  Using  Wavelet 
Packet  Transforms 


In  this  chapter,  we  approach  the  problem  of  plant  diagonalization  by  seeking 
wavelet  packet  transforms  which  approximate  the  SVD  factors  of  the  plant  matrix. 
We  consider  systems  whose  plant  matrix  G  is  a  fixed,  real  valued  matrix.  We  also 
assume  that  G  is  square  and  has  full  rank.  Such  a  model  is  suitable  for  “quasi¬ 
static”  linear  systems,  i.e.  stable  linear  systems  whose  transient  response  can  be 
neglected. 

The  idea  of  using  wavelet  transforms  in  the  context  of  feedback  for  large  control 
networks  comes  from  the  work  of  Chou,  Guthart,  and  Flamm  [14].  These  authors 
have  shown  that  for  a  special  class  of  linear  systems,  the  wavelet  transform  can  be 
used  to  transform  the  input  and  output  vectors  into  bases  which  allow  for  the  effi¬ 
cient  computation  of  a  feedback  control  law.  Further,  the  multiresolution  property 
of  the  wavelet  transform  leads  to  a  natural  implementation  on  a  hierarchical  con¬ 
trol  network.  Motivated  by  this  work,  we  seek  to  use  transforms  from  the  wavelet 
packet  as  the  basis  transformations  to  achieve  plant  diagonalization. 
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We  begin  the  discussion  with  the  naive  suggestion  that  the  matrix  factorization 
technique  of  singular  value  decomposition  can  be  used  for  plant  diagonalization. 
The  SVD  factors  accomplish  the  task  of  diagonalization,  but,  as  we  will  show, 
they  are  not  feasible  to  implement  on  a  large  control  network  because  they  are 
computationally  intensive  and  centralized. 

We  then  move  on  to  introduce  some  of  the  basic  concepts  of  wavelet  theory 
with  an  eye  toward  using  wavelet  transforms  for  plant  diagonalization.  We  con¬ 
centrate  on  the  implementation  of  transforms  from  the  wavelet  packet  and  show 
that  these  transforms  are  computationally  efficient  and  have  a  natural  parallel 
implementation  on  a  control  network. 

Next,  we  show  how  to  search  the  wavelet  packet  to  find  the  transforms  which 
most  closely  diagonalize  the  plant  matrix.  These  transforms  are  approximations 
of  the  SVD  factors.  An  algorithm  to  search  the  wavelet  packet  is  borrowed  from 
Wickerhauser  [43].  Motivated  by  this  work,  we  develop  a  cost  function  suitable  for 
the  purpose  of  matrix  diagonalization.  This  cost  function  is  then  used  to  search 
the  wavelet  packet  for  the  SVD  approximations. 

Finally,  we  demonstrate  this  idea  with  two  examples.  In  the  first  example, 
wavelet  packet  transforms  are  used  to  approximately  diagonalize  the  plant  matrix 
of  a  system  composed  of  a  flexible  membrane  driven  by  a  linear  array  of  electrostatic 
actuators.  The  second  example  looks  at  the  approximate  diagonalization  of  the 
plant  matrix  of  a  flexible  beam  driven  at  resonance. 

2.1  Singular  Value  Decomposition 

We  wish  to  find  basis  transformations  which  diagonalize  the  real  valued  plant 
matrix,  G.  This  can  be  accomplished  using  the  well  known  matrix  factorization 
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technique  of  singular  value  decomposition  (SVD).  Here,  we  briefly  discuss  SVD 
and  demonstrate  how  it  can  be  used  in  the  context  of  our  problem.  A  complete 
description  can  be  found  in  Strang  [35]. 

Consider  the  real  valued  n  x  n  matrix  G.  SVD  states  that  G  can  be  factored 
such  that 

G  =  Q1HQT,  (2.1) 

where 

1.  The  columns  of  Qi  e  Hnx?l  are  the  eigenvectors  of  GGT . 

2.  The  columns  of  Q2  G  Rnxn  are  the  eigenvectors  of  GTG. 

3.  The  n  x  n  matrix  S  is  diagonal,  and  values  on  the  diagonal  are  the  square 
roots  of  the  eigenvalues  of  GGT  and  GTG. 

The  matrices  Q\  and  Q2  are  orthogonal,  i.e.  QiQj  =  Q2Q2  =  where  In 
is  the  n  x  n  identity  matrix.  The  set  containing  all  n  x  n  orthogonal  matrices  is 
denoted  0(n).  Multiplying  G  by  Qj  on  the  left  and  Q2  on  the  right  yields 

Q[GQ2  =  S.  (2.2) 

In  the  present  context,  G  is  a  plant  matrix  of  an  //-input,  n-output  system.  In 
other  words, 

V  =  Gu,  (2.3) 

where  u  is  the  //-dimensional  input  vector  and  y  is  the  //-dimensional  output 
vector. 

Consider  the  change  of  variables  y  =  Qjy  and  u  =  Q^u.  Substitution  into  the 
plant  equation  (Equation  1.1)  yields 

V  =  Q1PQ2U 
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(2.4) 


=  Em. 

Hence,  the  zth  element  of  the  transformed  input  vector  u  affects  only  the  ith 
element  of  the  transformed  output  vector  y  for  i  —  1,2, ...  ,n.  As  a  result,  the 
entire  system  can  be  controlled  using  local  feedback,  i.e.  each  input  is  determined 
based  solely  on  the  value  of  its  corresponding  output.  In  other  words,  the  SVD 
factors  achieve  the  stated  goal  of  plant  diagonalization. 

Unfortunately,  the  transforms  generated  by  SVD  are  computationally  intensive 
to  implement;  for  an  n  dimensional  output  vector  y,  it  takes  0(n2)  operations  to 
compute  y.  Additionally,  the  value  of  every  element  of  y  is  required  to  compute 
each  element  of  y.  This  fact  renders  the  decentralized  nature  of  the  resulting  basis 
meaningless.  Fortunately,  we  can  overcome  these  implementation  issues  by  using 
wavelet  packet  transforms  to  approximate  the  SVD  factors.  This  is  the  subject  of 
the  remaining  part  of  this  chapter.  In  the  next  section,  we  lay  the  groundwork  for 
this  idea  by  reviewing  some  of  the  principles  of  wavelets  and  the  wavelet  packet. 

2.2  Wavelets 

Wavelets  have  been  the  subject  of  much  work  in  recent  years  [36]  [15].  A  wavelet 
is  a  localized  function  which,  together  with  its  dilations  and  translations,  forms  a 
basis  for  a  specified  function  space.  In  other  words,  if  0(f)  is  a  suitable  “mother” 
wavelet  for  a  given  space,  then  any  function  /(f)  belonging  to  that  space  can  be 
written  as 

fit)  =  ^2  -  k)  (2.5) 

The  associated  wavelet  transform  is  the  map  which  takes  a  function  from  the 
original  space  and  returns  the  coefficients  of  the  wavelet  expansion,  a,jk,j,k  G  Z. 
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The  wavelet  transform  simply  gives  us  another  way  to  represent  the  information 
contained  in  the  function  f(t). 

Wavelets  have  a  large  number  of  characteristics  and  applications.  Below  we 
focus  on  a  few  that  we  will  use.  Since  we  will  be  using  the  wavelet  transform  on 
discrete  signals  of  finite  length  (vectors),  the  following  discussion  is  limited  to  that 
case. 


2.2.1  The  Fast  Wavelet  Transform 


One  of  the  most  important  features  of  wavelets  is  that  the  wavelet  transform  is 
extremely  easy  and  fast  to  implement  for  discrete  spaces.  In  fact,  for  discrete 
signals,  the  fast  wavelet  transform  (FWT)  can  is  implemented  as  an  iterated  bank 
of  two  low  order  digital  filters  combined  with  a  down  sampling  step.  The  FWT 
is  briefly  described  here.  A  complete  discussion  of  wavelets  and  the  connections 
between  wavelets  and  filter  banks  can  be  found  in  Strang  and  Nguyen  [36]. 

The  Haar  wavelet  transform  can  be  implemented  as  an  iterated  bank  of  digital 
filters  as  shown  in  Figure  2.1.  This  filter  bank  is  called  the  “analysis  bank”.  The 
the  lowpass  filter  C  and  the  highpass  filter  D  are  written  in  ^-transform  form  as 


C(z) 

D(z) 


1  1  -i 

V2  +  W 

s/i  Vi 


(2.6) 

(2.7) 


The  output  of  each  filter  is  fed  into  a  downsampling  step,  denoted  in  Figure  2.1 
as  l  2.  The  output  of  the  downsampler  is  the  input  vector  with  every  other  element 
removed. 

The  output  of  the  highpass  channel  of  the  filter  bank  is  kept;  these  numbers 
are  part  of  the  wavelet  transform.  The  output  of  the  lowpass  channel  is  fed  into 
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Figure  2.1:  Filter  bank  implementation  of  the  fast  wavelet  transform. 

another  identical  filter  bank  and  the  process  is  iterated  until  the  desired  resolution 
is  reached.  In  the  case  of  a  vector  of  dimension  2n,  n  iterations  of  the  filter 
bank  completely  transform  the  vector;  the  output  of  each  rightmost  highpass  and 
lowpass  channel  is  a  single  real  number  so  that  no  further  filtering  is  possible. 

2.2.2  The  Wavelet  Packet 

The  wavelet  packet  is  a  collection  of  orthogonal  transforms  which  can  be  imple¬ 
mented  using  filter  banks.  Where  the  fast  wavelet  transform  is  realized  by  repeat¬ 
edly  passing  the  output  of  the  lowpass  filter  into  another  filter  bank,  the  entire 
wavelet  packet  is  realized  by  repeatedly  passing  the  outputs  of  both  channels  to 
the  next  filter  bank.  This  filtering  process  is  represented  in  Figure  2.2.  The  packet 
formed  using  the  filter  bank  from  the  Haar  wavelet  transform  is  called  the  Haar- 
Walsh  Packet.  For  a  vector  of  length  n  (an  even  power  of  two),  this  filtering  process 
is  iterated  log2(n)  times  so  that  the  output  of  each  filter  at  the  end  of  the  tree  (the 
right  of  the  tree  as  depicted  in  Figure  2.2)  is  a  single  real  valued  number. 

This  tree  of  filters  contains  a  large  number  of  orthonormal  representations  of 
the  input  vector.  For  example,  the  coefficients  taken  from  the  n  filters  at  the  end 
of  the  tree  provide  a  representation  of  the  input  in  a  basis  analogous  to  the  Short 
Time  Fourier  Transform  [36].  Additional  bases  can  be  found  by  “pruning”  the  tree, 
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Figure  2.2:  Filter  bank  implementation  of  the  wavelet  packet. 

i.e.  removing  filter  banks  from  the  tree.  The  rule  for  this  pruning  is  simple:  a  filter 
bank  block  can  only  be  removed  if  all  banks  which  use  the  output  of  that  bank  have 
already  been  removed.  Associated  with  each  basis  is  an  orthogonal  basis  transform 
which  is  realized  by  the  resulting  configuration  of  filter  banks.  The  wavelet  packet 
is  the  collection  of  all  such  transforms. 

Another  useful  representation  of  a  wavelet  packet  for  a  finite  length  vector  is 
shown  in  Figure  2.3.  Notice  that  each  level  of  the  tree  has  twice  as  many  blocks 
as  the  level  before.  The  blocks,  however,  are  each  half  as  long,  leaving  the  same 
number  of  coefficients  in  each  level.  If  we  define  a  “graph”  of  blocks  to  be  any 
set  of  blocks  in  the  tree  such  that  any  vertical  line  passes  through  exactly  one 
block,  then  the  collection  of  coefficients  from  any  graph  gives  a  representation  of 
the  input  vector  with  respect  to  some  orthonormal  basis  [43]. 

The  wavelet  packet  provides  a  large  library  of  orthonormal  basis  transforms. 
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level  L 


Figure  2.3:  Tree  representation  of  wavelet  packet. 

The  entire  packet  can  be  computed  in  0(nlog(n))  operations,  so  it  is  feasible  to 
create  such  a  library  even  for  large  n.  Further,  the  resulting  library  is  arranged 
in  a  tree  form.  As  we  will  see,  this  allows  the  library  to  be  searched  using  fast 
recursive  algorithms. 

2.2.3  Inverse  Transforms 

The  FWT  has  an  inverse  which  can  also  be  implemented  as  an  iterated  bank  of 
filters,  as  shown  in  Figure  2.4.  Since  this  filter  bank  reconstructs  the  original  signal 
it  is  called  the  “synthesis  bank” .  Here,  the  inputs  are  upsampled  before  they  are 
passed  into  the  filters.  Upsampling,  denoted  by  |  2,  inserts  a  zero  between  every 
element  of  the  input  vector.  For  the  case  of  the  Haar  wavelet,  the  filters  C*  and 
D*  are  written  in  ^-transform  form  as 

1  1  -i 

— 7=  +  —j=Z 

V2  V2 


C*(z) 
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Figure  2.4:  Filter  bank  implementation  of  the  inverse  fast  wavelet  transform. 


Figure  2.5:  Filter  bank  implementation  of  a  wavelet  packet  transform  along  with 
its  inverse. 

1  1  -i 

~j=  +  ~j=Z 

y/2  y/2 

A  similar  synthesis  bank  can  be  used  to  invert  any  transform  from  the  wavelet 
packet.  In  this  case,  the  synthesis  filter  bank  must  be  a  mirror  image  of  the 
analysis  filter  bank,  with  (7,  D ,  and  j.  2  replaced  by  C *,  D*,  and  t  2,  respectively. 
An  example  of  a  wavelet  transform  and  its  inverse  in  filter  bank  implementation 
is  shown  in  Figure  2.5. 


D*(z) 
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2.2.4  Higher  Order  Filter  Banks 

The  previous  three  sections  show  the  filter  bank  implementation  of  the  FWT, 
wavelet  packet  transform  (WPT),  and  their  inverses  for  the  Haar  wavelet.  There 
are  many  other  wavelets  which  can  be  implemented  using  the  same  filter  bank 
structure  with  wisely  chosen  FIR  filters  C,  D,  C*  and  D*.  Clearly,  the  filters 
must  be  chosen  so  that  the  synthesis  bank  actually  inverts  the  analysis  bank.  This 
feature  is  sometimes  referred  to  as  “perfect  reconstruction”.  Other  conditions, 
such  as  orthogonality,  can  also  be  imposed  on  the  filters  to  yield  good  wavelet 
properties. 

An  important  set  of  suitable  filters  has  been  introduced  by  Daubechies  [15].  In 
addition  to  perfect  reconstruction,  Daubechies’  filters  are  orthogonal  and  satisfy 
the  “maximum  flatness”  condition,  which  means  that  the  frequency  responses  of 
the  filters  are  as  flat  as  possible  at  u  =  0  and  u>  =  tt.  Other  choices  of  filters  which 
yield  good  wavelets  are  discussed  by  Strang  and  Nguyen  [36]. 

2.2.5  FWT  Implementation  on  a  Control  Network 

The  FWT  or  any  other  transform  from  the  wavelet  packet  can  be  implemented 
very  efficiently  on  a  distributed  control  network.  The  key  to  this  is  to  exploit  the 
recursive  nature  of  these  transforms.  At  each  filter  bank  iteration,  only  local  in¬ 
formation  is  used.  This  local  information  is  processed  (filtered)  and  the  important 
information  is  passed  on.  This  process  is  repeated  until  the  transform  is  complete. 
We  take  advantage  the  locality  of  the  information  at  any  step  in  the  transform  to 
distribute  the  necessary  calculations  over  a  control  network  connected  to  form  a 
hierarchy.  We  first  discuss  the  implementation  of  the  Haar- Walsh  wavelet  packet. 
Then  we  discuss  a  possible  extension  to  include  more  general  wavelet  packets. 
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These  ideas  will  be  made  more  specific  in  an  example  in  Section  2.4. 


The  Haar  Wavelet 


Here,  we  show  how  to  use  a  hierarchical  control  network  can  be  used  to  compute 
the  FWT  for  a  vector  of  length  8  using  the  Haar  wavelet.  Once  the  transform  is 
understood  for  an  8-vector,  extension  to  other  vector  lengths  is  trivial.  Extensions 
to  other  bases  in  the  Haar- Walsh  packet  will  be  briefly  discussed  at  the  end  of  this 
section  and  illustrated  by  the  example  in  Section  2.4. 

As  we  have  already  discussed,  the  FWT  using  the  Haar  wavelet  can  be  imple¬ 
mented  using  the  filter  bank  in  Figure  2.1  where  the  filters  C  and  D  are  given  by 
Equations  2.6  and  2.7.  After  downsampling,  the  output  of  the  lowpass  part  of  the 
first  filter  bank  is 

x(2)  +  x(l) 
x(4)  +  x(3) 
x(6)  +  x(5) 
x(8)  +  x(7) 

The  output  the  first  highpass  filter  after  down  sampling  is 


(2.10) 


Cx  =  —= 

V2 


Dx  = 


1 

71 


x(2)  —  x(l) 
x(4)  —  x(3) 
x(6)  —  x(5) 
x(8)  —  x(7) 


(2.11) 


Note  that  each  element  of  Cx  and  Dx  rely  on  adjacent  elements  of  x.  The  same 
is  true  if  we  look  at  the  outputs  of  the  next  filter  bank: 


C{Cx) 


J_  (Cx){2)  +  (Cx)(l) 
72  (Cx)(4)  +  (Cx)(3) 


(2.12) 
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and 


D(Cx )  =  -= 

V  '  V2 

Likewise,  the  outputs  of  the  final  filter  bank  are 

C(CCx)  =  — =  [{CCx){  2)  +  (CCx){  1)]  (2.14) 

v  2 

and 

D(C'CV)  =  -^=  [(CCx)(2)  -  (CCx)(l)} .  (2.15) 

v  2 

So  the  computation  of  C  and  D  rely  on  adjacent  elements  of  the  previous  filter 
bank’s  output  at  every  level.  This  is  the  locality  we  will  take  advantage  of  with  a 
hierarchy. 

Consider  the  tree  depicted  in  Figure  2.6.  The  nodes  at  the  bottom  level  contain 
the  elements  of  the  input  vector  x.  This  information  is  passed  up  to  the  next  level 
where  the  individual  elements  of  C  and  D  are  calculated  separately.  The  elements 
of  Cx  are  then  passed  to  the  next  level.  The  process  is  then  repeated  until  we 
reach  the  top  of  the  tree.  The  transform  is  then  given  by  the  coefficients  in  Dx, 
DCx ,  DCCx ,  and  CCC(x). 

The  inverse  FWT  can  be  implemented  on  a  hierarchy  by  using  the  same  idea, 
only  in  reverse.  Here,  we  start  at  the  top  of  the  tree  and  work  down,  using  the 
filters  C*  and  D*  instead  of  C  and  D. 

Implementing  the  FWT  and  its  inverse  on  this  hierarchy  has  two  big  advan¬ 
tages.  First,  communication  is  kept  to  a  minimum  since  each  node  receives  only 
the  information  it  needs  to  compute  its  part  of  the  transform.  Also,  the  connection 
scheme  required  is  simple:  each  node  needs  to  communicate  with  at  most  three 
other  nodes  (its  parent  and  its  two  children).  Second,  a  large  number  of  the  neces¬ 
sary  computations  are  done  in  parallel,  at  the  same  time,  on  different  processors. 


(CV>( 2)  -  (Cz)(l) 

(Cx)( 4)  -  (Cx)( 3) 


(2.13) 
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Figure  2.6:  Tree  implementation  of  the  Haar  wavelet  transform. 

For  example,  it  takes  a  total  of  12  operations  to  compute  Cx  (8  multiplications,  4 
additions).  Assuming  each  of  these  calculations  takes  one  cycle  to  complete,  this 
means  it  would  take  a  single  processor  12  cycles  to  compute  Cx.  Consider  instead 
the  hierarchy  approach.  It  takes  only  3  operations  to  compute  1  element  of  Cx. 
Since  each  element  of  Cx  is  calculated  on  its  own  processor,  all  of  the  elements 
can  be  computed  simultaneously.  Hence,  in  the  hierarchy  implementation,  it  takes 
only  3  cycles  to  compute  all  of  Cx. 

For  these  reasons,  we  propose  to  implement  the  FWT  on  a  control  network 
connected  to  form  a  hierarchy.  Recall  that  each  network  node  has  its  own  processor, 
so  massive  parallel  processing  capability  exists.  The  data  from  a  given  node’s 
onboard  sensor,  at  any  point  in  time,  can  be  thought  of  as  one  element  of  the 
input  vector,  x.  Our  goal  is  to  create  a  control  network  which  can  emulate  the 
hierarchy  discussed  above. 


23 


Figure  2.7:  Connection  scheme  to  realize  hierarchy  implementation  of  the  Haar 
wavelet  transform  using  a  distributed  control  network. 

One  such  possibility  is  shown  in  Figures  2.7  and  2.8.  Figure  2.7  depicts  a 
connection  scheme  for  a  linear  array  of  eight  nodes.  This  array  conld  be,  for 
example,  an  array  of  intelligent  strain  sensors  mounted  along  a  flexible  beam. 
Figure  2.8  shows  how  this  connection  scheme  can  be  used  to  implement  the  FWT. 
This  “hierarchy”  uses  a  limited  amount  of  communication  and  takes  advantage  of 
its  parallel  processing  abilities  in  computing  the  FWT. 

A  similar  process  can  be  used  to  implement  any  other  basis  transformation  in 
the  Haar- Walsh  packet.  In  this  case,  the  computations  are  made  as  in  the  wavelet 
transform  case:  each  node  can  compute  its  elements  of  the  filtered  vector  based 
on  information  from  its  two  children.  The  difference  lies  in  the  communication. 
Instead  of  just  passing  the  output  of  the  lowpass  filter  up  to  the  next  level,  it  may 
be  necessary  to  pass  the  output  of  the  highpass  filter  or  the  outputs  of  both  filters. 

Higher  Order  Wavelets 

In  the  previous  section  we  showed  how  to  efficiently  implement  any  transform  in  the 
Haar- Walsh  packet  on  a  hierarchical  control  network.  Since  the  Haar  filters  look  at 
two  adjacent  elements  of  the  input  vector  at  a  time  and  the  filter  is  followed  by  a 
down  sampling  of  two,  a  simple  dyadic  tree  was  the  natural  choice  for  the  network 
topology.  For  more  general  wavelets,  such  as  the  family  of  wavelets  proposed 
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Figure  2.8:  Hierarchy  implementation  of  Haar  wavelet  transform  using  a  dis¬ 
tributed  control  network. 
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Figure  2.9:  Multihierarchy  required  to  realize  hierarchy  implementation  of  a  dis¬ 
crete  wavelet  transform  with  4th  order  filter. 

by  Daubechies  [15],  the  filters  look  at  n  adjacent  elements  of  the  input  vector. 
The  resulting  output  is  then  downsampled  by  a  factor  of  two.  This  leads  to  a 
“multihierarchy”  topology.  Figure  2.9  shows  the  structure  of  such  a  hierarchy  for 
a  wavelet  filter  with  4  coefficients. 


2.3  Approximation  of  SVD  Factors  with  Wavelet 
Packet  Transforms 

In  the  previous  two  sections  we  have  laid  the  groundwork  for  the  main  result  of 
this  chapter.  In  Section  2.1  we  showed  that  the  SVD  factors  of  the  plant  matrix 
can  be  used  to  diagonalize  the  plant,  but  they  are  not  feasible  to  implement  on  a 
large  control  network.  In  Section  2.2  we  showed  that  wavelet  packet  transforms 
could  be  efficiently  implemented  on  a  control  network,  but  we  gave  no  indication 
of  how  they  could  be  used  to  diagonalize  a  plant  matrix.  Here,  we  combine  these 
two  ideas.  The  approach  is  straightforward.  The  wavelet  packet  contains  a  large 
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number  of  orthogonal  transforms.  Each  SVD  factor  is  also  orthogonal.  To  End  an 
“approximate”  SVD  transform,  we  simply  search  the  wavelet  packet  to  End  the 
transform  which  is,  in  some  sense,  closest  to  the  SVD  factor. 

This  idea  is  borrowed  from  Wickerhauser  [43],  who  proposed  this  approach  to 
perform  principal  component  analysis  on  a  set  of  random  data  vectors.  Using 
the  fact  that  the  Karhunen-Loeve  basis  is  the  minimum  entropy  basis  of  a  given 
data  set,  Wickerhauser  was  able  to  devise  a  cost  function  which  allowed  for  an 
efficient,  recursive  search  of  the  wavelet  packet.  A  similar  approach  can  be  used 
to  approximate  the  SVD  factors  of  the  plant  matrix. 

2.3.1  The  Cost  Function 

We  assume  that  the  plant  matrix  G  G  !Rn x n  is  square  and  full  rank.  Finding 
a  suitable  SVD  factor  Q[  is  equivalent  to  finding  the  0  G  O(n)  such  that  0T 
diagonalizes  the  symmetric  matrix  GGT .  Towards  this  goal,  we  develop  a  cost 
function  Ji(0)  which  is  globally  minimized  when  the  matrix  QGGTQT  is  diagonal. 
We  begin  by  writing 

£=[<&<&•••,<£],  (2.16) 

where  g{,  g%, . . . ,  are  the  columns  of  G.  Now  define 

xe  =  e  g  =  [©<,;,  e*  •  ■  ■ ,  e9‘] .  (2.17) 

Let 

Me  =  QGGTQT  =  XeXl.  (2.18) 

The  fill  diagonal  element  of  M©  is  then  given  by 

n 

[*fe]«  =  J>e]„[xa, 

3= 1 
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(2.19) 


n 

= 

3= 1 

Since  both  G  and  0  are  nonsingnlar,  each  [Mq]u  is  strictly  positive.  We  now  define 
a  cost  function  on  0(n)  by  taking  the  product  of  the  diagonal  elements  of  M©: 

n 

Me)  t  UlMel 

i=  1 
n  n 

=  nE(e»J)?-  (2-20) 

*= 1  j=l 

We  justify  that  Ji(0)  is  a  suitable  cost  function  to  be  minimized  with  the  following 
claim: 

Claim  2.3.1  Let  the  matrix  Qi  G  0{n )  be  a  left  SVD  factor  of  the  full  rank  matrix 
G  G  Rnxn?  i.e.  Qi  diagonalizes  GGT .  Then 

MQi)  <  Ji(0)  (2.21) 

for  all  0  G  0(n).  Furthermore,  Ji(Qj)  =  Ji(0)  if  and  only  if  0T  diagonalizes 
GGT. 


In  order  to  prove  this  claim,  we  require  some  matrix  facts.  We  state  these  facts 
without  proof  [35]. 


Fact  2.3.1 


det 


/ 

A  B 

\ 

V 

C  D 

/ 

det(A)det(D  -C A"1  B), 


(2.22) 


where  A  is  a  square,  invertible  matrix,  D  is  a  square  matrix,  B  and  C  are  matrices 
with  dimensions  compatible  to  A  and  D. 


Fact  2.3.2  A  symmetric  A  is  positive  definite  if  and  only  if  every  square  sub¬ 
matrix  in  the  upper  left  corner  of  A  has  positive  determinant. 
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Now  we  are  ready  to  prove  the  claim. 

Proof  of  Claim: 

First,  we  derive  an  expression  for  J\ ( Qj ) •  Note  that 


mqt 


Qt1GGtQ1 

a\  0  •  •  •  0 

0  a\  •••  0 


(2.23) 


where  cr©  <T2,.  . . ,  on  are  the  singular  values  of  G.  Since  G  is  full  rank,  all  of  its 
singular  values  are  strictly  positive.  Now 


MQl)  =  n[%] 

i— 1 
n 

=  ip? 

i= 1 

=  det  (GGt).  (2.24) 


Note  that  det(GG'T)  =  det(0GGrT0T)  =  det (M©)  for  any  0  G  0{n).  Hence,  in 
order  to  show  that  </i(Qf)  <  Ji(0)  it  suffices  to  show  that 


n 


Dehne  Tn  =  M©. 
7fc+i  G  R  such  that 


det(M©)  <  n^]«.  (2.25) 

i= 1 

Given  rfc+1  G  Rfc+lxfe+1,  dehne  Tk  eHkxk,  bk+1  G  Rfc  and 


Tfc+i  — 


(2.26) 


|_  bk+ 1  7fc+i  J 

for  k  =  n  —  l,n  —  2,..., 2,1.  Note  that  Tx  G  R.  Let  7 1  =  Note  that 
7^  =  [M<s]kk  and  each  >  0.  Since  M©  >  0,  Fact  2.3.2  tells  us  that  >  0  for 
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each  k.  Fact  2.3.1  then  gives  the  following  expression  for  det(Ffc): 


det(rfc)  =  det(rfc_i)det(7fc  -  blYk\bk) 

=  det(rfc_1)(7fc-6^r^16fc),  (2.27) 

k  =  2,3 , ,n.  We  can  use  this  equation  to  recursively  find  the  following  expres¬ 
sion  for  det  (Mq): 

det(M©)  =  det(rn) 

=  det(rn_i)(7„  —  b^T^\bn) 

=  det(rn_2)(7n_i  -  ^_ir^26„_i)(7„  -  bn) 


det(ri)n(7fc  -blTk\bk) 


k= 2 


=  7i 


(2.28) 


k= 2 


Note  that  each  bkTk\bk  >  0  since  Tfc_i  >  0  =7  F^  >  0.  Using  Equation  2.27 
and  the  fact  that  det(rfc)  >  0  for  each  k ,  we  see  that  (7*.  —  bkTk\bk)  >  0  for  each 
k  —  2,  3, . . . ,  n.  This  gives  the  inequality 


n  n 

71  -  hkTk-ibk )  <  n*- 


(2.29) 

k= 2  k= 1 

The  positive  dehniteness  of  each  guarantees  that  this  inequality  will  be  strict 
unless  all  of  the  off-diagonal  elements  of  M©  are  0.  Conversely,  if  all  of  the  off- 
diagonal  elements  of  M©  are  zero,  the  the  expression  in  Equation  2.29  becomes  an 
equality.  H 

Hence,  we  can  End  a  suitable  Qj  by  searching  0(n )  for  the  global  minimum  of 
./, .  It  is  convenient  to  redefine  J\  as 


■m©)  =  £  log  ( E  (0^)  ■ 


(2.30) 


i= 1 
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In  a  similar  manner,  we  can  define  a  cost  function  JyfB)  which  can  be  used  to 
search  for  QT2. 

n  /  n  \ 

^(0)=E'og  E(e»;),2  h  <2-31) 

i=l  \i=l  / 

where  </[,  <?£, . . .  grn  are  the  transposes  of  the  rows  of  G. 

2.3.2  Searching  the  Wavelet  Packet 

In  this  section  we  address  the  problem  of  searching  the  wavelet  packet  to  find  the 
“best”  approximations  of  the  SVD  factors  Q\  and  Q%.  We  find  these  approxima¬ 
tions  by  minimizing  the  functions  </i(@)  and  J2(0)  for  ©  £  L.  where  L  denotes  the 
collection  of  orthogonal  transforms  contained  in  the  wavelet  packet.  The  0  €  L 
which  minimizes  J\  gives  us  the  approximation  Q  j  while  the  0  which  minimizes 
J2  yields  the  approximation  Q\.  ()[  is  the  wavelet  packet  transform  which  most 
nearly  diagonalizes  the  matrix  GGT  and  Q?>  is  the  wavelet  packet  transform  which 
most  nearly  diagonalizes  the  matrix  GTG.  The  algorithm  used  to  search  the  packet 
is  due  to  Wickerhauser  [43]. 

The  transform  Q\  is  the  solution  of  the  minimization  problem 

min  </].(©),  (2.32) 

een 

where 

n  /  n  \ 

J1(°)  =  El°g  E(00.2  •  (2  33) 

i= 1  \j= 1  / 

To  understand  the  recursive  search  developed  by  Wickerhauser,  we  refer  back  to 
the  tree  representation  of  the  wavelet  packet  shown  in  Figure  2.3.  Any  transform 
in  the  wavelet  packet  can  be  represented  by  a  collection  of  blocks  on  the  tree  with 
the  property  that  any  vertical  line  drawn  through  the  tree  intersects  exactly  one 
block.  Such  a  collection  is  called  a  graph.  Let  Q  be  a  given  graph  and  let  0  be 
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the  associated  wavelet  packet  transform.  If  the  elements  of  a  vector  g  are  placed 
in  the  top  block  (block  1,  level  1),  then  the  transformed  vector  Qg  is  given  by 
the  coefficients  located  in  the  blocks  of  the  graph  Q.  The  blocks  in  any  graph 
are  independent  of  each  other  in  the  following  sense:  for  any  given  block  in  the 
graph,  the  other  blocks  in  the  graph  can  be  changed  to  create  a  new  graph  without 
affecting  the  coefficients  contained  in  the  original  block. 

Note  that  the  ith  component  in  the  outer  summation  in  the  right  hand  side 
of  Equation  2.33  depends  only  on  the  ith  element  of  each  transformed  column 
of  the  plant  matrix,  Qg?,  j  =  1,2, ...  ,n.  This  allows  us  to  construct  a  tree  of 
costs.  We  begin  by  constructing  a  tree  representation  for  each  g?,  j  =  1,2,...,  n. 
We  then  square  each  element  in  each  of  these  trees.  Next,  we  add  the  n  squared 
trees  together  on  an  element  by  element  basis  to  create  one  tree  called  the  “sum 
of  squares”  tree.  Finally,  the  sum  of  squares  tree  is  converted  to  the  cost  tree  by 
taking  the  logarithm  and  adding  the  resulting  numbers  within  each  block  together. 

The  cost  tree  assigns  a  cost  to  each  block  in  the  tree  representation  of  the 
wavelet  packet.  Let  0  be  a  given  wavelet  packet  transform  and  let  Q  be  the 
associated  graph.  The  cost  of  0  can  be  computed  by  adding  together  the  costs 
associated  with  the  blocks  contained  in  Q.  The  cost  computed  in  this  manner  is 
identical  to  the  cost  described  by  Equation  2.33. 

Now  the  cost  tree  can  be  efficiently  searched  to  End  the  wavelet  packet  trans¬ 
form  which  yields  the  lowest  cost.  This  is  done  by  comparing  the  cost  of  each  block 
with  the  sum  of  the  “best  costs”  of  its  two  child  blocks.  Of  course,  finding  the  best 
cost  of  each  child  block  requires  comparing  the  child  to  the  child’s  children,  and 
so  on.  This  recursion  terminates,  and  the  collection  of  blocks  that  results  gives 
the  wavelet  packet  transform  that  yields  the  lowest  cost.  The  resulting  transform 
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is  then  used  as  the  SVD  factor  approximation  Q\ .  The  complexity  of  this  search, 
including  the  construction  of  the  cost  tree,  is  0(jr\ogp),  where  p  is  the  number 
of  rows  in  the  plant  matrix.  A  more  rigorous  discussion  of  the  properties  of  this 
search  can  be  found  in  [43]. 

A  similar  search  using  J2(@)  is  used  to  obtain  the  SVD  factor  approximation 

Ql 

2.4  Examples 

Here  we  introduce  two  example  systems  to  demonstrate  approximate  diagonaliza- 
tion  using  wavelet  packets.  The  first  system  is  a  simplified  model  of  a  flexible 
membrane  driven  by  a  linear  array  of  point  forces.  This  example  is  motivated  by 
the  MEMS  mirror  array  developed  by  Bifano  et.al.  [3].  The  second  example  is  that 
of  a  flexible  cantilever  beam  driven  at  resonance.  These  examples  will  also  be  used 
to  demonstrate  the  effectiveness  of  other  ideas  which  will  be  presented  later  in  the 
thesis. 

2.4.1  Simplified  Flexible  Membrane 

The  system  we  consider  has  32  inputs  and  32  outputs.  The  inputs  correspond 
to  electrostatic  actuators,  which  we  assume  act  as  point  forces  on  the  membrane. 
The  actuators  are  evenly  spaced  over  the  length  of  the  array.  Let  this  interval 
be  denoted  by  A.  The  outputs  correspond  to  the  displacement  of  the  mirror 
measured  at  the  positions  of  the  32  actuators.  We  assume  that  the  displacement 
Hi  due  to  the  input  Ui  is  simply  Ui.  To  determine  the  displacement  i/j  due  to  Ui,  we 
assume  that  the  deflection  of  the  mirror  at  the  positions  A  to  the  left  of  the  first 
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Figure  2.10:  Response  of  the  membrane  to  the  ith  input. 


actuator  and  A  to  the  right  of  the  32nd  actuator  is  always  zero.  The  displacement 
changes  linearly  between  the  two  end  points  and  the  position  of  the  ith  actuator 
(see  Figure  2.10).  We  also  assume  that  the  displacement  at  a  given  position  is  the 
snm  of  the  displacements  due  to  each  of  the  actuators  so  that  the  plant  can  be 
represented  by  a  matrix.  The  plant  matrix  G  is  depicted  in  Figure  2.11.  Here,  the 
intensities  of  the  pixels  correspond  to  the  values  of  the  entries  in  the  matrix. 

We  applied  the  method  described  in  this  chapter  to  find  the  transforms  from 
the  Haar-Walsh  packet  which  most  closely  diagonalize  the  plant  matrix  G.  The 
resulting  “approximately  diagonal”  matrix  is  shown  in  Figure  2.12.  The  diagonal- 
ization  of  the  plant  matrix  using  SVD  is  shown  in  Figure  2.13.  Comparing  the 
two,  we  see  that  the  approximately  diagonalized  matrix  resembles  the  matrix  di¬ 
agonalized  with  SVD.  A  more  quantitative  discussion  of  these  results  can  be  found 
in  Chapter  6. 

This  example  can  be  used  to  further  illustrate  the  implementation  of  the  wavelet 
packet  transforms.  In  the  tree  representation,  the  wavelet  packet  transform  which 
corresponds  to  the  output  transformation  Q]  is  shown  in  Figure  2.14.  In  this 
figure,  the  elements  of  the  transformed  data  vector  are  contained  in  the  solid  black 
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Plant  Matrix  (1 D  Membrane) 


Figure  2.11:  Plant  matrix  for  flexible  membrane. 


Approximately  Diagonalized  Plant  Matrix 


Figure  2.12:  Approximately  diagonalized  membrane  plant  matrix  using  a  transform 
from  the  Haar- Walsh  packet. 
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Diagonalized  Plant  Matrix  (Using  SVD) 


Figure  2.13:  Diagonalized  membrane  plant  matrix  using  exact  SVD  factors. 

boxes.  A  network  hierarchy  which  could  be  used  to  implement  this  transform  is 
shown  in  Figure  2.15.  The  circles  represent  numbers,  with  the  circles  in  the  leftmost 
column  representing  the  elements  of  y.  The  solid  circles  represent  elements  of  the 
transformed  output  vector.  A  pair  of  circles  with  a  box  around  them  represents 
a  computation  step.  The  node  on  which  the  computation  step  is  performed  is 
indicated  by  the  vertical  position  of  the  box;  if  the  box  is  in  the  same  row  with  yi: 
then  the  corresponding  computation  step  is  executed  on  the  ith  node.  In  such  a 
computation  step,  the  numbers  X\  and  x2  are  transformed  to  yield  (xi  +  x2) / y/2 
and  (— x\  +x2)/^/2.  The  lines  with  arrows  on  them  indicate  the  communication  of 
one  number  from  one  node  to  another.  The  transform  proceeds  from  left  to  right 
in  the  diagram. 

The  input  transformation,  which  in  this  example  happens  to  be  the  inverse  of 
the  output  transformation,  is  implemented  by  reversing  this  process.  To  visualize 
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Figure  2.14:  Tree  representation  of  the  wavelet  packet  basis  which  approximately 
diagonalizes  the  plant  matrix. 

how  the  inverse  transform  can  be  implemented  on  a  distributed  control  network  we 
again  look  to  Figure  2.15.  In  this  case,  the  numbers  we  start  with  are  the  elements 
of  u.  The  job  of  the  transform  is  to  transform  u  to  the  actual  input  vector  u.  The 
elements  of  u  are  indicated  in  the  figure  by  the  solid  circles.  At  the  beginning  of  the 
transform  they  are  distributed  among  various  nodes  on  the  network.  The  transform 
executes  across  the  figure  from  right  to  left  with  the  data  transfers  flowing  against 
the  directions  of  the  arrows.  As  in  the  output  transformation  figure,  the  boxes 
represent  a  computation  step  performed  on  the  specified  node.  In  this  case,  the 
numbers  x\  and  x2  are  transformed  to  yield  Jrx2)/\f7i  and  (x\  —  *2)/V2.  After 
this  computation  is  complete,  the  resulting  numbers  are  then  sent  to  the  next 
two  blocks  of  nodes  on  the  hierarchy.  The  sums  are  sent  to  the  upper  blocks  and 
the  differences  are  sent  to  the  lower  blocks.  When  the  transform  is  complete,  the 
elements  of  u  are  distributed  on  the  network  nodes  with  Ui  residing  on  the  ith 
node. 
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Figure  2.15:  Graphical  depiction  of  the  implementation  of  the  wavelet  packet  trans¬ 
form  shown  in  Figure  2.14. 
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2.4.2  Flexible  Beam  at  Resonance 


In  this  example  we  consider  a  flexible  cantilever  beam  driven  by  PZT  actuators 
as  described  in  Appendix  B.  Most  of  the  energy  in  the  response  of  the  beam 
lies  within  narrow  bands  in  the  frequency  domain  which  correspond  to  the  reso¬ 
nant  frequencies  of  the  beam.  Because  the  beam  naturally  rejects  non-resonant 
disturbances  very  well,  the  resonant  frequencies  are  separated  by  relatively  large 
frequency  bands  where  the  gain  is  very  small.  In  this  sense,  the  resonances  are 
isolated  from  one  another.  In  most  applications  where  feedback  control  is  used  for 
resonant  systems,  it  is  only  necessary  to  apply  contol  in  the  resonant  frequency 
bands. 

One  method  of  rejecting  near-resonant  disturbances  in  a  SISO  flexible  beam 
is  the  so-called  positive  position  feedback  (PPF)  method  developed  by  Fanson 
and  Caughey  [18].  This  technique  uses  a  bank  of  bandpass  filters  to  seperate  the 
output  into  its  resonances.  Each  resonance  is  then  controlled  independently.  A 
similar  “resonance  by  resonance”  approach  can  be  taken  for  a  MIMO  beam.  At 
resonance  (or  at  any  fixed  frequency)  a  MIMO  beam  can  be  modeled  by  a  fixed 
complex  matrix,  i.e. 

y  =  Gcu.  (2.34) 

The  matrix  Gc  is  just  the  transfer  function  matrix  of  the  beam  evaluated  at  the 
resonant  frequency.  The  complex  matrix  Gc  can  be  diagonalized  using  the  same 
basis  transforms  which  diagonalize  the  corresponding  magnitude  matrix  G.  i.e.  the 
matrix  containing  the  magnitudes  of  the  elements  of  Gc.  The  complex  elements 
of  the  vectors  y  and  u  represent  the  magnitude  and  phase  of  a  sinusiod  at  the 
resonant  frequency. 

Here  we  consider  a  flexible  beam  with  32  inputs  and  32  outputs  at  its  6th 
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Flexible  Beam  Plant  Matrix  (6th  Resonance) 
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Figure  2.16:  Plant  matrix  for  flexible  beam  at  6th  resonance. 

resonant  frequency.  The  magnitude  plant  matrix  for  this  system  is  depicted  in 
Figure  2.16.  The  wavelet  packet  based  approximate  diagonalization  of  the  magni¬ 
tude  matrix  is  depicted  in  Figure  2.17.  The  diagonalization  generated  using  exact 
SVD  factors  is  shown  in  Figure  2.18. 
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Chapter  3 


Recursive  Orthogonal  Transforms 


In  this  chapter  we  introduce  a  class  of  linear  operators  that  are  well  suited  to 
being  implemented  in  a  distributed  fashion  on  a  control  network.  The  operators  we 
present  have  two  features  which  make  them  particularly  useful  for  implementation 
on  a  control  network.  First,  the  required  computations  are  naturally  distributed 
to  take  advantage  of  the  parallel  processing  capability  of  the  network.  Second,  the 
transforms  can  be  chosen  so  that  only  nearest  neighbor  communication  is  required, 
a  fact  which  is  important  for  systems  with  limited  communication  bandwidth. 

We  begin  our  discussion  with  a  review  of  the  necessary  mathematical  tools  and 
concepts  that  will  be  required  here  and  in  following  chapters.  We  then  move  on 
to  an  ad  hoc  discussion  about  implementing  data  transformations  on  the  input 
and  output  vectors  of  one  and  two  dimensional  arrays  of  sensor /actuator  pairs. 
We  then  generalize  these  ideas  and  define  the  main  tool  which  will  be  used  in  the 
remainder  of  the  thesis:  recursive  orthogonal  transforms  (ROTs). 
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3.1  Preliminaries 


Here  we  present  the  mathematical  tools  and  concepts  necessary  to  develop  the 
discussion.  We  first  review  some  basic  properties  of  Lie  groups  and  vector  fields  on 
Lie  groups.  Next  we  turn  to  linear  algebra  for  a  brief  discussion  of  singular  value 
decomposition  and  symmetric  matrix  diagonalization.  Finally,  we  present  some 
miscellaneous  matrix  facts  which  will  be  required  later  in  this  document. 

3.1.1  Vector  Fields  on  Lie  Groups 

Here  we  present  some  requisite  information  on  Lie  groups  and  their  associated 
vector  fields.  A  comprehensive  treatment  of  this  subject  can  be  found  in  the 
references  [40],  [20],  [42],  and  [31].  A  nice  overview  appears  in  the  appendix  of 
[22],  We  concentrate  on  the  orthogonal  group,  0(n),  the  special  orthogonal  group, 
SO(n),  and  the  unitary  group,  U(n). 

Smooth  Manifolds 

Formally,  a  smooth  manifold  is  a  Hausdorff  topological  space  which  is  locally  dif- 
feomorphic  to  Rn  around  each  point.  An  example  of  a  smooth  manifold  is  the 
circle,  S'1  =  {(xi,x2)  G  R2  :  x\  +  x\  =  1}.  Around  each  point,  a  neighborhood  of 
the  circle  can  be  “flattened  out”  using  a  smooth,  invertible  mapping.  Hence,  the 
circle  is  locally  diffeomorphic  to  R.  With  this  intuition,  we  state  a  more  rigorous 
definition  of  a  smooth  manifold. 

Definition  3.1.1  A  smooth  manifold  is  a  topological  Hausdorff  space  M  along 
with  a  collection  of  charts  A  =  {(4>a,  Ua,  Va)  :  a  G  A}.  For  each  a  G  A,  Ua  and 
Va  are  open  subsets  of  M  and  Rn,  respectively,  and  </>«  :  Ua  — >  Va  is  a  smooth, 
invertible  map.  In  addition,  the  charts  must  satisfy  the  following: 
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1.  M  —  UaeA  Ua,  be.  the  charts  cover  M. 

2.  4>p  o  0”1  :  <j>a(Ua  D  Up)  C  Va  — >  (pfj{Ua  D  Up)  C  Vp  is  a  diffeomorphism  for 
each  a,  (3  G  A,  i.e.  the  charts  are  compatible  where  their  domains  intersect. 

3.  Any  chart  (0,  U,  V )  which  is  compatible  with  every  chart  in  A  is  also  con¬ 
tained  in  A,  i.e.  A  is  maximal. 

The  collection  A  is  called  an  atlas.  The  integer  n  is  called  the  dimension  of  the 
manifold.  The  collection  of  charts  allows  us  to  represent  an  open  neighborhood  of 
any  point  on  the  manifold  as  an  open  subset  of  ]Rn.  This  representation  is  referred 
to  as  the  local  coordinate  representation  of  M. 

Lie  Groups 

Here  we  introduce  Lie  groups,  which  are  a  special  type  of  smooth  manifold. 

Definition  3.1.2  A  group  is  a  set  G  along  with  an  operation  o  which  has  the 
following  properties: 

1.  g  o  h  G  G  for  all  g,  h  G  G,  i.e.  the  group  is  closed  under  the  operation  o. 

2.  (,9i  o  g2)  o  g3  =  g1  o  (g2  o  g3)  for  all  91,92, 93  e  G,  i.e.  the  operation  o  is 
associative. 

3.  There  exists  e  G  G  such  that  goe  =  eog  =  g  for  all  g  G  G,  i.e.  the  group 
contains  an  identity  element. 

4-  For  each  g  G  G,  there  exists  g~l  G  G  such  that  g  o  g~l  =  e,  i.e.  each  group 
element  has  an  inverse  which  is  also  contained  in  the  group. 

If  H  is  a  subset  of  G  which  is  also  a  group,  then  H  is  called  a  subgroup  of  G. 
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To  simplify  the  notation,  we  will  drop  the  use  of  the  o  symbol  so  that  go  h  will  be 
written  gh. 

Definition  3.1.3  A  Lie  group  G  is  a  smooth  manifold  with  a  group  structure  such 
that  the  group  operation  o  :  G  x  G  — >  G,  (g,h)  (->■  gh  is  a  smooth  map.  If  H  is  a 
subgroup  of  G  it  is  called  a  Lie  subgroup. 

We  are  most  interested  in  the  matrix  Lie  groups  0(n),  SO(n),  and  U(n).  Here 
we  define  these  groups  and  state  some  of  their  basic  properties  without  proof. 

The  Orthogonal  Group,  0(n): 

Abstractly,  the  orthogonal  group  can  be  thought  of  as  the  set  of  rotations  and 
reflections  in  Hn.  This  group  can  be  represented  as  the  following  matrix  Lie 
group: 

0{n)  =  {0e  ]Rnxn|©T©  =  1}  .  (3.1) 

O(n)  is  a  smooth  manifold  with  dimension  n(n  —  l)/2. 

The  Special  Orthogonal  Group,  SO(n): 

It  is  readily  seen  that  each  element  of  0(n )  has  determinant  equal  to  either  +1  or 
—  1.  The  subset  of  matrices  in  0(n)  with  determinant  equal  to  +1  is  also  a  Lie 
group.  This  subgroup  is  called  the  special  orthogonal  group  and  is  denoted  SO(n). 
The  special  orthogonal  group  can  be  thought  as  the  set  of  rotations  in  ]Rn.  SO(n ) 
has  the  same  dimension  as  O(n). 

The  Unitary  Group,  U{n): 

The  unitary  group  is  the  generalization  of  0(n )  to  complex  valued  matrices.  Specif¬ 
ically, 

U{n)  =  {0G  (Cnxn|0H0  =  1}  .  (3.2) 

The  real  dimension  of  U(n)  is  n2.  U{n )  contains  both  0(n )  and  SO[n). 
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Tangent  Spaces  and  Vector  Fields 


At  each  point  on  a  smooth  manifold  there  is  an  associated  vector  space  called  the 
tangent  space.  Intuitively,  if  we  imagine  an  infinitely  small  object  which  moves 
around  on  the  manifold,  the  tangent  space  at  a  given  point  contains  the  set  of  all 
possible  velocities  the  object  can  attain  at  that  point.  The  tangent  space  of  M  at 
the  point  x  is  denoted  TXM.  In  the  example  of  the  circle,  the  tangent  space  at  a 
given  point  x  is  the  line  tangent  to  the  circle  at  x,  or 

TxS1  =  {yeU2\yTx  =  0}.  (3.3) 


To  define  more  rigorously  what  we  mean  by  tangent  vectors  and  tangent  spaces, 
we  begin  by  defining  the  notion  of  a  curve  on  the  manifold.  Specifically,  a  curve 
through  x  G  M  is  a  smooth  function  /  :  (— e,  e)  — »  M,  some  real  e  >  0,  with 
/( 0)  =  x.  Let  (0,  U,  V )  be  a  coordinate  chart  with  x  G  U.  Let  Fx  denote  the  set  of 
all  possible  curves  at  x.  In  a  neighborhood  of  x,  the  curve  /  can  be  represented  as 
a  curve  in  ]Rn  using  local  coordinates,  4>o  f  :  (—5,  h)  — *  V  C  1R”  for  some  S  >  0.  If 
we  think  of  the  argument  of  /  as  time,  then  the  velocity  of  a  point  moving  along 
the  curve  at  the  point  x  is  given  in  local  coordinates  by  derivative  of  <po  f  at  x  =  0. 
Define  the  set 


=>/(*»)}. 


(3,4) 


[f]x  is  the  equivalence  class  of  all  curves  through  x  which  have  the  same  velocity 
at  x.  This  leads  to  the  following  definition: 


Definition  3.1.4  A  tangent  vector  to  M  at  the  point  x  is  an  equivalence  class 
[f }x,  where  f  G  Fx.  The  tangent  space  of  M  at  x,  TXM ,  is  the  set  of  all  tangent 
vectors  to  M  at  x.  The  tangent  bundle  of  M,  denoted  TM,  is  the  set  containing 
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the  tangent  spaces  of  every  x  E  M , 


TM  =  {TxM\x  E  M}  .  (3.5) 

An  element  of  the  tangent  bundle  is  an  ordered  pair  ( x,p )  where  x  is  a  point  on 
the  manifold  and  p  is  a  vector  in  TXM.  There  is  a  natural  projection  from  the 
tangent  bundle  to  the  manifold,  7 t  :  TM  — »  M;  (x,p)  H >  x.  We  are  now  ready  to 
define  a  vector  held: 

Definition  3.1.5  A  vector  held  on  a  smooth  manifold  M  is  a  map  X  :  M  -E  TM 
with  the  condition  that  n  o  X  is  the  identity  map. 

Hence,  the  vector  held  assigns  a  tangent  vector  px  to  each  point  on  the  manifold. 

The  flow  of  a  vector  held  starting  from  the  point  Xo  E  M  is  the  curve  x(-)  : 
1R  — »  M  that  satishes  the  differential  equation 

x(t)  =  px,  x(0)  =  x0.  (3.6) 


The  Lie  Bracket  and  Lie  Algebras 

Every  Lie  group  G  has  associated  with  it  a  vector  space  called  a  Lie  algebra  denoted 

0- 

Definition  3.1.6  A  Lie  algebra  is  a  vector  space  V  together  with  the  Lie  bracket 
[•,•]:  V  x  V  — ?•  V  such  that 

1.  [a,  b]  =  —[b,  a]  for  all  a,b  G  V . 

2.  [aa  +  fib,  c]  =  a  [a,  c]  +  fi[b,  c]  for  all  a,  fi  E  R,  a,b,c  E  V. 

3.  [a,  [b,  c]]  +  [c,  [a,  6]]  +  [b,  [c,  a]]  =  0  for  all  a,b,c  E  V  (The  Jacobi  Identity). 
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The  Lie  algebra  of  a  Lie  group  G  is  given  by  the  tangent  space  at  the  identity, 
0  =  TeG,  together  with  an  appropriate  Lie  bracket.  The  Lie  algebras  of  0(n), 
SO(n ),  and  U(n)  are  denoted  o(n),  so(n),  and  u(n),  respectively.  They  are 

o(n)  =  so(n)  =  {fle  !Rnxn|ft  =  -0T}  (3.7) 

and 

u (n)  =  {0  G  (Dnxn|0  =  -nH}  .  (3.8) 

The  Lie  bracket  for  o  (n),  so  (n),  and  u(n)  is  given  by  the  matrix  commutator, 

[A,  B}  =  AB-  BA.  (3.9) 

One  of  the  nice  properties  of  Lie  groups  is  that  the  tangent  space  at  the  identity 
can  be  translated  to  the  tangent  space  at  g  G  G  via  left  multiplication  by  g.  Hence, 
for  matrix  Lie  groups,  the  tangent  space  at  g  G  G  can  be  written 

TgG  =  {gn\n  e  g}  .  (3.10) 

Each  vector  0  in  the  matrix  Lie  algebra  naturally  generates  a  vector  held  on  the 
Lie  group  via  the  definition  Xq,  :  G  — *  TG ;  g  (->•  ( g ,  gQ).  Such  a  vector  held  is  said 
to  be  left  invariant.  Flows  of  Xq  are  curves  which  satisfy 

g  =  gtt,  g(0)  =  g0e  G.  (3.11) 

On  a  matrix  Lie  group  this  equation  is  a  linear  time  invariant  matrix  ODE.  The 
solution  is  given  by 

g(t)  =  g0eQt  (3.12) 

where  the  matrix  exponential  is  dehned  by  the  series 

Ak 

k= 0 
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Riemannian  Metrics,  Directional  Derivatives,  and  Gradient  Flows 

Each  smooth  function  /  :  M  — >  1R  has  associated  with  it  a  gradient  ascent  vector 
held.  At  each  point  in  M,  the  gradient  vector  held  assigns  a  tangent  vector  which 
points  in  the  direction  of  “steepest  ascent”,  i.e.  the  direction  in  which  /  increases 
the  most.  Following  [22],  we  present  the  tools  necessary  to  formalize  this  notion 
in  the  setting  of  smooth  manifolds. 

Definition  3.1.7  A  Riemannian  metric  on  a  manifold  M  is  a  collection  of  nonde¬ 
generate  inner  products  {(•,  -)x  :  TXM  — *  R}  for  each  x  G  M  where  (•,  -)x  depends 
smoothly  on  x.  A  manifold  with  such  a  specified  collection  is  called  a  Riemannian 
manifold. 

For  0{n)  and  SO(n),  the  tangent  space  at  any  given  point  is  a  subset  of  Rnx?l.  As 
a  result,  the  standard  inner  product  on  n  x  n  real  matrices,  (A,  B)  =  tr  (ArI?), 
provides  a  suitable  Riemannian  metric: 

(v)  :TgO(n)  xTgO(n)  ->  R 

(gtti,gSh)  ^  tr(f l^gTgQ2) 

=  tr(nffl2),  (3-14) 

where  G  o(n).  Similarly,  U(n)  inherits  a  Riemannian  metric  from(Dnxn: 

(y)  :  TgU{n)  x  TgU{n)  ->  R 

(gili,gil2)  ^  retr  (ftf gH g^) 

=  retr  (Q^Q2)  ,  (3.15) 

where  fli,  G  u(n)  and  retr  (•)  denotes  the  real  part  of  the  trace. 
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Consider  a  smooth  function  /  :  M  — >  R.  The  directional  derivative  of  /  at  a 
point  x  G  M  is  a  map 

Dfx  :  TXM  ->  R.  (3.16) 

We  will  use  Dfx{£)  to  denote  the  directional  derivative  of  /  at  x  in  the  direction 
£  e  TXM. 

Now  we  consider  the  gradient  vector  held  of  /,  which  is  denoted  as  V/.  The 
gradient  vector  held  is  dehned  as  the  vector  held  which  satishes  the  following  two 
properties: 

1.  V/(x)  G  TXM  for  all  x  G  M. 

2.  £>/*(£)  =  <V/(xU>  for  all  £  G  TXM. 

The  hrst  property  states  that  the  gradient  vector  is  contained  in  the  tangent  space 
of  M  at  every  point  in  M.  The  second  property  is  a  compatibility  condition  be¬ 
tween  the  gradient  and  the  directional  derivative.  For  a  given  Riemannian  metric, 
the  gradient  vector  held  which  has  both  properties  is  unique.  The  ascent  direction 
gradient  how  from  initial  condition  xq  G  M  is  given  by  the  solution  to 

x  =  V/(x),  x(0)  =  xo ■  (3-17) 

In  this  document  we  will  be  interested  in  gradient  hows  which  evolve  on  the 
compact  Lie  groups  SO{n )  and  U(n).  Here  we  state  some  useful  properties  of  a 
gradient  how  on  a  compact  manifold  without  proof  [22], 

Proposition  3.1.1  Consider  the  gradient  flow  given  by  Equation  3.17  where  the 
manifold  M  is  compact.  Then  for  any  x0  G  M : 

1.  The  solution  x(t)  exists  for  all  time  t  G  R. 
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2.  The  function  f(x(t ))  is  a  nondecreasing  function  oft. 


3.  The  solution  x(t)  converges  to  x*  G  M  such  that  V  f(x*)  =  0. 

3.1.2  Matrix  Facts 

Here  we  present  some  facts  about  matrices  which  will  be  required  for  upcoming 
calculations. 

Fact  3.1.1  Any  A  G  (Dnxn  can  be  written  as 

A  =  S  +  Tl ,  (3.18) 

where  S  is  hermitian  ( S  =  SH )  and  (12)  is  skew-hermitian  (12  =  —ClH).  Specifi¬ 
cally,  S  and  12  are  given  by 

S  =  ]-(A  +  AH)  (3.19) 

Zu 

12  =  ^  (A  -  Ah)  .  (3.20) 

This  fact  can  be  easily  verified  by  adding  the  expressions  for  S  and  12.  Clearly, 
this  fact  can  be  applied  to  A  G  ]Rnxn  by  replacing  the  hermitian  transpose  with 
the  regular  matrix  transpose. 

Fact  3.1.2  Let  S  =  SH  G  (Dnxn  and  12  =  -12 H  e  (Dnxn.  Then  retr(SQ)  =  0. 

Proof: 

We  begin  by  computing 

retr  ((ST2)H)  =  retr  (12 H SH) 

=  retr  (—US') 

=  —retr  (SCI) .  (3-21) 
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Since  retr  (A)  =  retr  (AH)  for  any  complex  matrix  A,  we  then  have  retr  (Sfi)  = 
— retr  (Sfi).  This  can  only  be  true  if  retr  (Sfi)  =  0.  gg 

Clearly,  if  S'  is  a  real  valued  symmetric  matrix  and  fi  is  real  valued  and  skew- 
symmetric,  then  tr  (Sfi)  =  0. 

Fact  3.1.3  Let  A  e  (Dnxn  be  written  A  =  S  +  fl  where  S  =  SH  and  fl  =  —flH. 
Let  A  =  —  AH  e(Dnxn.  Then  retr  (Ah)  =  retr(fi  A). 

This  fact  follows  directly  from  Fact  3.1.1  and  Fact  3.1.2. 

3.2  Parallel  Signal  Processing  on  Control  Net¬ 
works 

Here  we  provide  an  ad  hoc  discussion  about  the  implementation  of  coordinate 
transforms  on  the  input  and  output  vectors  of  a  large  control  network.  This  dis¬ 
cussion  leads  us  to  a  collection  of  orthogonal  transforms  which  have  a  special 
structure  which  allows  the  transforms  to  be  efficiently  implemented  on  the  sys¬ 
tems  we  consider.  We  restrict  the  discussion  to  control  networks  whose  nodes  are 
sensor /actuator  pairs.  We  first  investigate  the  case  where  the  nodes  are  distributed 
spatially  along  a  linear  array.  We  then  address  a  similar  control  network  whose 
nodes  are  distributed  on  a  rectangular  array. 

3.2.1  Implementation  on  1— D  Array 

In  this  section  we  address  the  problem  of  implementing  input  and  output  data 
transforms  on  a  control  network  where  n  sensor /actuator  pairs  are  distributed 
spatially  to  form  a  one  dimensional  array.  The  sensor/actuator  pairs,  or  nodes, 
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are  indexed  from  left  to  right;  the  leftmost  node  has  the  index  1.  For  each  i,  the 
ith  node  contains  one  sensor,  whose  output  is  yt ,  and  one  actuator,  whose  input 
is  Ui .  We  assume  that  each  node  on  the  array  can  communicate  with  both  of  its 
neighbors.  At  the  endpoints  of  the  array,  it  is  assumed  that  the  first  node  can 
communicate  with  the  last  node.  For  i  —  1,  2, . . . ,  n  —  1  we  say  that  the  (i  +  l)th 
node  is  to  the  right  of  the  ith  node  and  the  ith  node  is  to  the  left  of  the  (i  +  l)th 
node.  Similarly,  we  say  that  the  first  node  is  to  the  right  of  the  last  node  and  the 
last  node  is  to  the  left  of  the  first  node.  The  presented  data  transformation  method 
requires  only  nearest  neighbor  communication  and  distributes  the  computations 
required  for  the  transforms  in  a  natural  way. 

First  we  concentrate  on  the  problem  of  transforming  output  vector  y.  The  result 
of  this  operation  is  the  transformed  output  vector  y  =  Qiy,  where  the  orthogonal 
transform  Q\  has  a  special  structure.  We  address  only  the  implementation  of 
Q i,  we  do  not  address  what  Q\  actually  does  to  the  data.  This  question  will 
be  addressed  in  Chapter  4.  We  then  show  how  a  similar  method  can  be  used  to 
transform  the  input  vector,  yielding  u  =  Q^u. 

Transforming  the  Output  Vector 

Here  we  present  a  data  transformation  which  can  be  easily  applied  to  the  output 
vector  of  a  control  network.  We  call  the  transformed  data  vector  y  and  the  data 
transformation  that  results  is  denoted  by  Qi,  so  y  —  Qiy.  The  sensors  are  dis¬ 
tributed  along  a  line  and  each  sensor  has  its  own  microprocessor.  There  is  also 
an  actuator  paired  with  each  sensor.  The  combination  of  sensor,  actuator,  and 
microprocessor  will  be  referred  to  as  a  node.  The  ith  output  r/j  is  known  only 
to  the  ith  node.  The  object  is  to  compute  y  =  Qiy  while  limiting  the  required 
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communication  and  distributing  the  required  computation. 

The  implementation  we  propose  consists  of  a  sequence  of  steps  or  levels.  We 
will  present  these  steps  here  without  motivating  our  choices.  Our  reasoning  should 
become  clear  at  the  end  of  this  section.  Also,  to  simplify  this  presentation  we 
assume  that  n  is  even.  The  extension  to  odd  n  will  be  addressed  when  we  generalize 
these  ideas  in  Section  3.3. 

In  the  first  level,  each  evenly  indexed  node  passes  the  value  of  its  output  to  the 
node  to  its  left.  As  a  result,  for  each  odd  i,  the  ith  node  contains  yt  and  yi+i-  The 
ith  node  then  uses  its  local  microprocessor  to  compute 


&i 

-  0(*+ 1)/2 

Vi 

^i+1 

Vi+ 1 

(3.22) 


where  the  matrix  Oh+1y2  G  SO( 2)  for  each  odd  i.  The  output  of  the  first  level  is 
the  intermediate  vector  a,  which  can  be  written 


a  =  0i  y, 


(3.23) 


where 


0i 


e\  o 
o  e\ 


o 

0 


(3.24) 


L0  •  0  0nl  2J 

When  the  first  level  is  complete,  the  ith  node,  i  odd,  contains  the  numbers  a, 
and  cti+ 1.  To  begin  the  second  level,  the  ith  node  passes  at  to  the  node  to  its  left 
and  a*+i  to  the  node  to  its  right.  Recall  that  nth  node  is  to  the  left  of  the  first 
node,  so  aq  gets  sent  to  the  nth  node.  With  the  exception  of  the  nth  node,  the  ith 
node  for  each  even  i  contains  a*  and  a,;_,  i .  The  nth  node  contains  an  and  (i  \ .  Then 
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each  ith  node,  i  =  2, 4, 6, . . . ,  n  —  2,  uses  its  local  processor  to  compute 


The  nth  node  computes 


h/  - 1 

bi 

bn—  1 
bn 


Oi/ 2 


tti+l 


(3.25) 


=  ^/2 


ai 


(3.26) 


The  matrix  d)y2  G  SO( 2)  for  each  i  =  2, 4,  6, . . . ,  n.  The  output  of  the  second  level 
is  the  vector  b ,  which  can  be  written 


1 

'  01 

0 

1 - 

o 

b2 

= 

0 

df 

0 

0 

0 

o2n/  2 . 

01 

0 

0 

= 

0 

df 

0 

0 

0 
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a  2 

0j3 
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0 
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0  0  •••  0  1 


1  0 


0  0 


a\ 

&2 


=  e2Pea 


(3.27) 


Substituting  Equation  3.23  into  this  expression  yields 

b  =  02  PeOiy. 


(3.28) 


The  matrix  Pe  is  the  permutation  matrix  which  implements  a  circular  left  shift  of 
the  vector  a. 
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When  the  second  level  is  finished  the  zth  node  contains  6.,_i  and  for  each 
even  i.  To  start  the  third  level,  the  zth  node  (i  even)  sends  to  the  node  on 
its  left  and  it  sends  bt  to  the  node  on  its  right.  Recall  that  the  first  node  is  to  the 
right  of  the  nth  node,  so  bn  gets  sent  to  the  first  node.  Like  the  first  and  second 
steps,  each  oddly  indexed  node  uses  its  local  processor  to  rotate  the  elements  of  b 
which  it  contains,  yielding  the  vector  c: 


=  QsPob.  (3.29) 

Substituting  Equation  3.28  into  this  expression  yields 

c  =  Q3P0Q2PeQiy.  (3.30) 

This  process  of  local  rotations  and  permutations  is  repeated  for  L  levels,  some 
L>  0,  and  y  is  defined  to  be  the  output  of  the  Lth  level, 

V  =  ©L-Pl  '  '  '  ©2-p2©l?/)  (3.31) 


56 


Figure  3.1:  Graphical  description  of  parallel  output  processing  on  a  1-D  array  of 
smart  sensors.  Data  flows  from  the  bottom  up  in  the  diagram. 


where,  for  i  =  2, 3, . . . ,  L, 


{PQ  for  i  odd 
Pe  for  i  even 


This  process  is  depicted  graphically  in  Figure  3.1. 
The  resulting  transform  Q i  is  of  the  form 


Q 1  —  ©L-Pl  '  '  '  ©2^2©l- 


(3.32) 


(3.33) 


Clearly,  any  orthogonal  transform  which  has  this  structure  can  be  implemented  in 
a  manner  which  naturally  distributes  the  necessary  computations  and  eliminates 
the  need  for  global  communication. 
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Transforming  the  Input  Vector 


In  the  previous  section  we  introduced  an  output  vector  transformation  Qi  that 
could  be  implemented  on  the  network  using  a  series  of  2  x  2  decoupled  rotations 
and  permutation  steps.  Here  we  consider  the  coordinate  transformation  for  the 
input  vector,  Q2.  We  start  by  assuming  that,  like  Q\.  Q2  is  of  the  form 

Q2  =  &lPl  • '  •  (3.34) 

where  the  dys  are  block  diagonal  rotation  matrices  of  the  same  form  as  the  (-).,s 
which  compose  Q\.  To  simplify  the  discussion,  we  assume  that  the  number  of 
levels,  in  the  transform  Q2  is  L,  the  number  of  levels  in  the  transform  Qi.  We  also 
assume  that  L  is  odd.  The  extension  to  more  general  cases  will  be  considered  in 
Section  3.3. 

In  order  to  understand  how  the  input  transformation  is  implemented,  we  first 
recall  that  the  proposed  controller  implementation  is  based  on  plant  diagonaliza- 
tion.  In  this  form,  the  controller  performs  three  basic  tasks:  it  transforms  the  data 
vector  y  to  y,  it  chooses  the  transformed  input  vector  u  based  on  y,  and  it  trans¬ 
forms  u  to  create  the  actual  input  vector  applied  to  the  plant,  u.  At  the  end  of  the 
output  transformation  step,  each  oddly  indexed  node  contains  two  elements  of  the 
transformed  output  vector  y.  In  particular,  the  ith  node  contains  the  transformed 
outputs  y,;  and  yp+i)-  1  As  seen  from  the  transformed  input  and  output  vectors, 
the  plant  is  diagonal.  The  controller  matrix  can  then  also  be  diagonal.  Hence,  each 
Ui  can  be  computed  based  solely  on  the  corresponding  yt.  This  means  that  the 
controller  can  perform  its  second  task,  computing  u,  in  a  decentralized  manner  by 
having  the  ith  node,  i  odd,  compute  ut  and  up+i)  on  its  local  processor.  So  at  the 

1Remember:  we  have  assumed  that  L  is  odd;  if  L  was  even  then  the  outputs  would  be 
contained  in  the  evenly  indexed  nodes  and  the  ith  node  would  contain  yu~i)  and  yi. 
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beginning  of  the  input  transformation  step,  the  transformed  inputs  Ui  and  wp+i) 
are  contained  in  the  ith  node  for  each  odd  i.  The  goal  of  the  input  transformation 
step  is  to  compute  u  =  Q^u  in  a  distributed  manner  on  the  network  so  that  the 
input  Ui  is  applied  to  the  actuator  on  the  ith  node. 

Now  we  are  ready  to  describe  the  process  used  to  transform  u.  The  input  we 
wish  to  compute  is 

u  =  **P?*T---PZ*Iu.  (3.35) 

Like  the  output  transformation,  the  input  transformation  is  computed  in  levels. 
In  the  first  level,  the  intermediate  vector  a  =  Qju  is  computed.  The  matrix  (l>] 
is  block  diagonal  so  a  can  be  computed  in  a  decentralized  manner  on  the  oddly 
indexed  nodes  of  the  network.  Next,  the  permutation  matrix  P[  is  realized  using 
nearest  neighbor  communication  where  each  odd  node  sends  one  element  of  a  to 
each  of  its  two  evenly  indexed  neighbors.  The  vector  b  =  can  then  be 

computed  in  a  decentralized  manner  on  the  evenly  indexed  nodes  of  the  network. 
This  process  is  repeated  until  finally  the  expression 

u  =  ^P?^---PZ$tlu  (3.36) 

is  realized.  Operating  under  the  assumption  that  L  is  odd,  the  elements  of  u  will 
be  contained  in  the  evenly  indexed  nodes,  with  the  ith  node  containing  W(j_ p  and 
Ui .  The  final  step  of  the  input  transformation  process  is  for  each  evenly  indexed 
ith  node  to  send  the  value  of  up_i)  to  the  node  on  its  left.  The  result  is  that  u 
has  been  computed  and  each  Ui  resides  on  the  ith  node,  ready  to  be  applied  as  the 
input  to  the  ith  actuator. 


59 


3.2.2  Implementation  on  2— D  Array 

In  this  section  we  show  how  to  design  specially  structured  orthogonal  transforms 
Qi  and  d)2  to  be  implemented  on  a  control  network  where  the  sensor/actuator 
pairs  are  distributed  spatially  to  form  a  two  dimensional  rectangular  array.  Let 
nr  be  equal  to  the  number  of  rows  and  nc  be  equal  to  the  number  of  columns  in 
the  array.  To  simplify  the  presentation,  we  assume  that  nc  and  nr  are  both  even. 
The  nodes  in  the  array  are  numbered  left  to  right,  top  to  bottom.  The  node  in 
the  ith  row  and  the  j th  column  is  numbered  j  +  nc(i  —  1).  Depending  on  which 
is  more  descriptive,  we  will  refer  to  this  node  as  either  the  (i,j) th  node  or  the 
(j  +  nc{i  —  l))th  node.  The  difference  should  be  clear  from  context. 

The  (i,  j)th  node  is  said  to  be  double  odd  if  both  i  and  j  are  odd.  Likewise,  the 
(i,  j)th  node  is  said  to  be  double  even  if  both  i  and  j  are  even. 

The  spatial  arrangement  of  the  array  naturally  divides  the  nodes  on  the  array 
into  three  classes:  interior  nodes,  edge  nodes,  and  corner  nodes.  As  the  name 
implies,  the  corner  nodes  are  at  locations  (1, 1),  (1  ,nc),  (nr,  1),  and  (nr,nc).  Nodes 
which  are  in  the  first  row,  last  row,  first  column,  or  last  column  but  are  not  corner 
nodes  are  referred  to  as  edge  nodes.  Specifically,  we  will  refer  to  these  edges  as  top 
edge,  bottom  edge,  left  edge,  and  right  edge,  respectively.  The  remaining  nodes 
are  interior  nodes. 

Array  Connections 

Here  we  assume  that  the  nodes  on  the  array  are  connected  to  form  a  torus.  To 
achieve  this  connection  topology,  imagine  that  the  top  and  bottom  edges  of  the 
array  are  brought  together  so  that  the  array  forms  a  cylinder  with  the  left  and 
right  edges  forming  the  circles  at  either  end  of  the  cylinder.  The  cylinder  is  then 
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stretched  and  bent  to  bring  the  left  and  right  edges  (now  circles)  together,  forming 
a  torus.  Remember:  this  torus  represents  the  topology  of  the  connections;  the 
physical  layout  of  the  array  has  not  changed. 

Each  node  on  the  array  has  eight  neighbors  on  the  torus.  For  interior  nodes,  the 
(z,  j)th  node  has  left  and  right  neighbors  (nodes  ( i ,  j  —  1)  and  (z,  j+ 1),  respectively), 
upper  and  lower  neighbors  (nodes  (z  —  1  ,j)  and  (z  +  1,  j),  respectively),  and  four 
diagonal  neighbors  (nodes  (z  —  l,j  —  1),  (z  — l,j  +  l),  (i  +  l,j  —  1),  and  (z  +  l,j  +  l). 
Because  of  toroidal  connection  topology,  nodes  on  the  edges  and  corners  have 
neighbors  on  the  other  side  of  the  array.  Modular  division  can  be  used  to  identify 
the  neighbors  of  a  general  element  on  the  array.  For  the  (i,j) th  node,  the  eight 
neighbors  are 


left  neighbor  ( z ,  ((j  —  2)  mod  nc )  +  1) 

right  neighbor  ( i ,  ((j)  mod  nc)  +  1) 
upper  neighbor  (((z  —  2)  mod  nr)  +  1  ,j) 

lower  neighbor  (((z)  mod  nr)  +  1,  j) 

upper-left  neighbor  (((z  —  2)  mod  nr)  +  1,  (( j  —  2)  mod  nc)  +  1) 

upper-right  neighbor  (((z  —  2)  mod  nr)  +  1,  ((j)  mod  nc)  +  1) 

lower-left  neighbor  (((z)  mod  nr)  +  1,  (( j  —  2)  mod  nc)  +  1) 

lower-right  neighbor  (((z)  mod  nr)  +  1,  (( j )  mod  nc)  +  1) 


(3.37) 


The  data  transformation  method  we  present  here  requires  each  node  to  have  a 
two-way  connection  with  a  subset  of  its  neighbors.  In  particular,  each  double  odd 
node  shares  a  connection  with  its  right  neighbor,  its  lower  neighbor,  and  each  of 
its  four  diagonal  neighbors.  This  connection  scheme  is  depicted  in  Figure  3.2. 
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•  • 


double  odd  nodes 
double  even  nodes 


Figure  3.2:  Required  connections  for  2-D  array. 


Transforming  the  Output  Vector 


Here  we  discuss  the  transformation  of  the  output  vector  y  on  a  2-D  array.  The 
resulting  transformed  output  vector  is  denoted  y  =  Qiy,  where  Q\  is  an  orthogonal 
transform  with  a  special  structure.  The  data  vector  y  has  dimension  nrnc.  Let 
y\i,j)  denote  the  output  of  the  sensor  on  the  (i,j)th  node.  The  output  of  the 
sensor  on  the  (i,j) th  node  is  contained  in  the  (j  +  nc(i  —  l))th  element  of  y,  i.e. 
y\i,j)  =  y(j+nc(i- 1)).  Conversely,  the  Mil  element  of  y  is 


Vk  = 


k_ 

nc 


+  1,  ((k  -  1)  mod  nc)  +  1  )  , 
where  |_mj  is  the  largest  integer  less  than  or  equal  to  x. 


(3.38) 


As  in  the  case  of  the  1-D  array,  the  implementation  we  propose  consists  of  a 
sequence  of  levels.  In  the  first  level,  each  double  odd  node  (i,j)  receives  the  values 
of  the  sensor  outputs  from  its  right  neighbor  (■ i,j  +  1),  left  neighbor  {i  +  1,  j),  and 


lower-right  neighbor  (i+  l,j  +  1).  As  a  result,  each  double  odd  node  (i,  j)  contains 
the  following  four  elements  of  the  output  vector,  y: 


=  y(j+nc(i- 1)) 

(3.39) 

y\i,j  +  !) 

=  y({j+l)+nc(i-l)) 

(3.40) 

y\i  +  i,i) 

=  y(j+nci) 

(3.41) 

y\i  +  i,j  +  i) 

=  y({j+l)+nci) 

(3.42) 

If  we  order  the  double  odd  nodes  from  left  to  right,  top  to  bottom,  then  the  (i,  j)th 
node  is  also  the  kth  double  odd  node,  where 


_  j  +  1  nci-  1 
~  2  +  2  2 


(3.43) 


Define  the  intermediate  data  vector  a  such  that  the  elements  of  y  contained  in  the 
kth  double  odd  node  are  equal  to  a^k- i)+i  through  a^k- 1)+4-  In  simple  language, 
this  means  that  the  first  four  elements  of  a  are  the  four  outputs  contained  in 
the  first  double  odd  node,  the  next  four  elements  are  the  four  outputs  contained 
in  the  second  double  odd  node,  and  so  on.  The  vector  a  is  a  permutation  of 
y.  Let  {ei,  e2,es, ,  enrHc}  denote  the  canonical  basis  vectors  for  Rn.  Then  we 
can  construct  the  permutation  matrix  I\  such  that  «  =  I\y  using  the  following 
algorithm: 


for  i  —  {1,  3, 5, . . . ,  nr  —  1} 
for  j  =  {1,3,5, . . .  ,nc  -  1} 
lot  h  =  i±l  -4-  acini 

lei  ft.  2^22 

let  the  (4 (k  —  1)  +  l)th  column  of  P1  be  eJ+nc(j_ i) 

let  the  (4 [k  —  1)  +  2)th  column  of  Pi  be  e(j+i)+ncp_i) 

let  the  (A(k  —  1)  +  3)th  column  of  Pi  be  e3+n<:l 
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let  the  (4 (k  —  1)  +  4)th  column  of  Pi  be  e(j+i)+?- 
end  j  loop 
end  i  loop 


This  communication  step  is  depicted  graphically  in  Figure  3.3. 

Next,  using  its  local  microprocessor,  each  double  odd  node  k  computes 


&4(fc-l)+l 

a4(fc-l)+l 

&4(fc-l)+2 

=  e\ 

a4(fc-l)+2 

&4(fc-l)+3 

«4(fc-l)+3 

&4(fc-l)+4 

a4(fc-l)+4 

(3.44) 


where  9l  G  S'0(4),  fe  =  1,2,3,...,  ncnr/A.  The  intermediate  vector  b  is  the  output 
of  the  first  level  and  can  be  written 


where 


©i  = 


b  =  QiPiy, 

e\  o 

0  e\  ••• 


(3.45) 


0 

0 


0 


®  ^ncTir/  4 


(3.46) 


In  the  second  level,  the  computations  will  take  place  on  the  double  even  nodes. 
We  number  the  double  even  nodes  left  to  right,  top  to  bottom  so  that  the  (i,j)th 
node  is  the  kth  double  even  node,  where 

j  ,  nci-  2 


k  = 


(3.47) 


2  2  2 

Note  that  the  kth  double  even  node  is  the  lower-right  neighbor  to  the  kth  double 
odd  node. 
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At  the  beginning  of  the  second  level,  each  double  even  node  receives  one  element 
of  b  from  each  of  its  four  diagonal  neighbors,  which  are  double  odd  nodes.  As 
depicted  in  Figure  3.3  this  communication  step  creates  the  new  intermediate  data 
vector  c.  Specifically,  we  set 


c4(fc-l)+l 

=  h(kul-l)+4 

(3.48) 

c4(fc-l)+2 

=  ^4(fcUr-l)+3 

(3.49) 

c4(fc-l)+3 

=  &4(fc„-l)+2 

(3.50) 

c4(fc-l)+4 

=  ^4(fc;r  — 1)+1  J 

(3.51) 

Where  kui ,  kur ,  ku,  and  kir  are  the  positions  of  the  double  odd  nodes  which  are, 
respectively,  the  upper-left,  upper-right,  lower-left,  and  lower-right  neighbors  of 
the  /cth  double  even  node.  Let  (i,j)  denote  the  coordinates  of  the  kth  double  even 
node.  Then 


kui 

u 

rvUr 

ku 

klr 


k 


j  mod  nc  +  2 
2 


+ 


nci  —  1 
~2  2 


j_ 

2 


nc  i  mod  nr 
~2  2 


j  mod  nc  +  2 

2 


+ 


nc  i  mod  nr 
~2  2 


(3.52) 

(3.53) 

(3.54) 

(3.55) 


The  intermediate  vector  c  is  a  permutation  of  b.  The  permutation  matrix  Pe  which 
satisfies  c  =  Peb  is  given  by  the  following  algorithm: 


for  i  —  { 2, 4,  6, ... ,  nr} 

for  j  =  {2,4,6,  ...,nc} 

let  fc  =  f  +  ir¥ 
let  k,,,i  =  k 


let  ktir  — 


_  j  mod  nc+ 2  nc  i- 2 


2  2 
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let 

hi  -- 

=  i  + 
2  ~ 

nc  i 

2 

mod  nr 
2 

let 

hr  " 

__  j  mod  r 

~  2 

ic+2  _|_  n 

',c  i  mod 

2  2 

nr 

let 

the 

m  - 

-i) 

+  l)th 

column 

of 

Pe 

be 

e4 (kuV 

— 1)+4 

let 

the 

m  - 

-i) 

+  2)th 

column 

of 

Pe 

be 

^  A (kur 

— 1)+3 

let 

the 

m  - 

-i) 

+  3)th 

column 

of 

Pe 

be 

eA(kn- 

-l)+2 

let 

the 

m  - 

-i) 

+  4)th 

column 

of 

Pe 

be 

emlr- 

-1)+1 

end  j  loop 
end  i  loop 


Next,  using  its  local  microprocessor,  each  double  even  node  k  computes 


^4(fe-l)+l 

c4(fc-l)+l 

<^4(fe-l)+2 

=  01 

c4(fc-l)+2 

<^4(fc-l)+3 

c4(fc-l)+3 

<^4(fc-l)+4 

c4(fc-l)+4 

(3.56) 


where  Q\  G  50(4)  for  each  k.  The  intermediate  vector  d  is  the  output  of  the  second 
level  and  can  be  written 


d=02Pee1P1y, 


(3.57) 


where 


e\  o  •  •  •  o 

o  el  ■  ■  •  o 


(3.58) 


0 


0  elcnr/A 


A  third  level  can  be  computed  in  a  similar  manner.  Here,  the  double  even  nodes 
pass  their  data  back  to  the  double  odd  nodes  which  perform  the  4D  rotations.  The 
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The  communication  patterns  for  the  first  3  levels  are  depicted  in  Figure  3.3. 
The  permutation  matrices  Pi,  Pe,  and  PQ  give  us  a  mathematical  representation 
of  these  communication  patterns. 

Levels  2  and  3  can  be  repeated  for  L  levels,  some  L  >  0,  and  y  is  defined  to  be 
the  output  of  the  Lth  level, 

y  =  ©l-Pl  •  •  •  ©2-P2G>i-Pl2/,  (3.60) 
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Level  1: 


Level  2: 


Level  3: 


•+o  m+o 

6\>  6^® 

9*0  9*0 

6\>  6^® 


Figure  3.3:  Communication  patterns  for  first  3  levels  of  parallel  output  transfor¬ 
mation  on  a  2-D  array  of  smart  sensors. 


where,  for  i  —  2, 3, . . . ,  L, 


{PQ  for  i  odd 
Pe  for  i  even 


Hence,  we  have  constructed  a  transform  of  the  form 


(3.61) 


Qi  =  ©lPl---02P20iPi 


(3.62) 


which  is  easily  implemented  in  a  distributed  manner  on  a  control  network. 


Transforming  the  Input  Vector 

Here  we  consider  the  problem  of  transforming  the  input  vector  on  a  2-D  array.  We 
begin  by  assuming  that  the  input  transformation  Q 2  is  of  the  same  form  as  Q\  in 
the  previous  section,  i.e. 


Q2  —  &lPl  ■  ■  ■  <L-P2<L-Pl, 


(3.63) 
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where  the  <3>jS  are  of  the  same  form  as  the  (-),;S  in  0 \ .  The  number  of  levels  in  Q 2 
is  equal  to  L,  which  is  the  same  as  the  number  of  levels  in  Q\ .  We  also  assume 
that  L  is  odd.  The  task  of  this  section  is  to  compute  the  corresponding  input 
transformation  u  =  Q^u  on  the  network.  As  in  the  case  of  the  1-D  array,  07  1 
can  be  implemented  on  the  2-D  array  by  reversing  the  steps  used  in  the  input 
transformation. 

To  begin,  recall  that  each  element  of  the  desired  transformed  input  vector  u 
resides  on  the  same  node  as  the  corresponding  element  of  the  transformed  output 
vector  y.  As  a  result,  u  is  distributed  over  the  network  with  each  double  odd  node 
containing  4  elements  of  u  (here  we  have  assumed  that  L  is  odd).  In  the  first  level, 
the  intermediate  data  vector  a  =  is  computed  on  the  double  odd  nodes.  Since 
is  a  block  diagonal  matrix,  a  can  be  computed  in  a  decentralized  manner.  Each 
double  odd  node  then  passes  one  element  of  a  to  each  of  its  diagonal  neighbors 
(the  double  even  nodes),  implementing 

b  =  Pea  (3.64) 

=  Pja  (3.65) 

=  (3.66) 

The  double  even  nodes  then  use  their  local  processors  to  compute 

c  =  (3.67) 

=  *Li  Pltfv-  (368) 

This  process  is  repeated  until  the  desired  input 

U  =  P-r*P^T...pT*T~  (3.69) 

is  achieved.  The  final  permutation  step  Pj  sends  the  elements  of  u  to  the  nodes 
on  which  they  will  be  applied. 
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3.3  Recursive  Orthogonal  Transforms 


In  the  previous  section  we  discussed  two  possible  methods  of  performing  data 
transformations  on  the  input  and  output  vector  of  an  array  of  smart  sensors  and 
actuators.  Here  we  provide  some  insight  into  our  choices  for  these  transforms. 
We  then  generalize  this  notion  to  define  what  we  call  the  recursive  orthogonal 
transform  (ROT). 

The  output  transforms  from  the  previous  section  are  all  of  the  form 

Q  =  QLPL...Q2P2Q1P1  (3.70) 


where,  for  i  =  1,2 , . . . ,  L,  the  OjS  are  block  diagonal  matrices  with  orthogonal 
blocks  on  the  diagonal  and  the  PjS  are  permutations  of  the  identity  matrix.  The 
represent  decoupled,  local  rotations  of  the  data  vector  performed  in  a  distributed 
fashion  on  the  control  network.  The  permutation  matrices  represent  a  “shuffling” 
of  the  data  vector  resulting  from  the  exchange  of  information  on  the  network. 

The  wavelet  packet  transforms  discussed  in  the  previous  chapter  can  be  written 
in  the  form  of  Equation  3.70.  To  see  this  connection,  consider  the  Haar  wavelet 
transform,  Q,  operating  on  a  4-dimensional  vector  y.  It  is  easy  to  verify  that 
the  wavelet  transform  can  be  written  as  the  following  matrix  sequence  of  matrix 
operations: 


Qy 


1 

1 

0 

0 

1 

1 

0 

0 

s/2 

s/2 

s/2 

s/2 

1 

1 

0 

0 

0 

0 

1 

1 

s/2 

s/2 

s/2 

V2 

0 

0 

1 

0 

1 

V2 

1 

s/2 

0 

0 

0 

0 

0 

1 

0 

0 

1 

s/2 

1 

s/2  J 

(3.71) 
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Employing  a  permutation  matrix,  this  expression  becomes 


Qy 


1 

s/2 

1 

s/2 

0 

0 

1 

0 

0 

0 

1 

s/2 

1 

s/2 

0 

0 

1 

vA 

1 

s/2 

0 

0 

0 

0 

1 

0 

1 

s/2 

1 

s/2 

0 

0 

0 

0 

1 

0 

0 

1 

0 

0 

0 

0 

1 

V2 

1 

s/2 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

1 

s/2 

1 

s/2  J 

(3.72) 


The  submatrix 

1  1 

V2  \/2 

_ J_ 

y/2  s/2 

is  orthogonal;  it  is  a  rotation  in  the  plane  by  -45  degrees.  Hence,  the  wavelet 
transform  is  a  special  case  of  a  transform  of  the  form  given  by  Equation  3.70.  The 
same  is  true  for  any  transform  in  the  wavelet  packet. 

Another  connection  to  the  this  idea  is  the  Euler  angle  representation  of  SO( 3), 
which  can  be  thought  of  as  the  group  of  orientations  of  a  rigid  body  in  3-dimen¬ 
sional  Euclidean  space  [30] .  The  Z YZ  Euler  angles  tell  us  that  any  orientation  can 
be  achieved  by  a  rotation  about  the  z  axis  by  an  angle  a,  followed  by  a  rotation 
about  the  new  y  axis  by  an  angle  (3,  followed  by  a  rotation  about  the  z  axis  by  an 
angle  7.  The  first  rotation  is  represented  by  the  matrix 


Rz(a) 


cos(a )  —sin(a)  0 
sin(a)  cos(a)  0 

0  0  1 


(3.74) 


This  matrix  is  block  diagonal,  and  the  blocks  on  the  diagonal  are  orthogonal.  The 
y  axis  rotation  can  be  written 


RM 


cos(/3 )  0  — sin((3 ) 
0  1  0 
sin(a)  0  cos(a ) 
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10  0  cos(P)  —sin(f3)  0  10  0 

0  0  1  sin(a)  cos(a)  0  0  0  1 

0  1  0  0  0  1  0  1  0 

=  P3RZ(P)P2.  (3.75) 

Now  any  0  G  SO( 3)  can  be  represented  by  the  composition  of  the  three  Euler 
rotations,  i.e. 

0  =  Rz{pi)Ry{p)Rz(a) 

=  Rz^)P3Rz(P)P2Rz(a)  (3.76) 

for  some  a,/3, 7  G  [0,27t).  Hence,  we  can  represent  any  0  G  SO( 3)  as  a  product  of 
block  diagonal  matrices  and  permutations  of  the  identity. 

Our  idea  is  to  extend  the  concept  of  Euler  angles  to  higher  dimensional  groups. 
This  leads  us  to  define  the  recursive  orthogonal  transform. 

Definition  3.3.1  A  recursive  orthogonal  transform  (ROT)  is  a  linear  operator  of 
the  form 

e  =  p1o1p2e2---pLeL,  (3.77) 

for  some  L,  where  for  each  i  —  1,  2, . . . ,  L, 

1.  The  matrix  0*  is  of  the  form 

e\  o  •  •  •  o 

0  9i  ■■■  0 

©<  =  ,  (3.78) 

.  0  •••  0  _ 

0j  G  SOiriij)  (or  0)  G  U (ni:j)), 

mi 

Uij  =  n,  (3.79) 

3= i 


72 


and  Os  represent  appropriately  dimensioned  blocks  of  zeros. 


2.  The  matrix  Pj  is  a  permutation  of  the  n  x  n  identity  matrix. 

The  integer  L  is  called  the  level  of  the  ROT.  The  structure  of  the  ROT  is  defined 
by  L,  Pi,  mi,  and  for  i  =  1,2, ...  ,L,  j  =  1,2, ... ,  mi.  These  quantities  are 
called  the  ROT  parameters.  The  ©j,  i  =  1, 2, . . . ,  L,  are  called  the  ROT  variables. 

In  the  remainder  of  this  thesis,  we  will  assume  that  the  ROT  parameters  are 
given  and  fixed.  We  will  address  the  problem  of  finding  the  ROT  variables  so  that 
the  resulting  ROT  most  closely  diagonalizes  a  given  plant  matrix. 
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Chapter  4 


Symmetric  Matrix 
Diagonalization  and  SVD  using 
ROT 


In  Chapter  1  we  discussed  the  possibility  of  applying  linear  data  transformations  to 
the  input  and  output  vectors  of  a  MIMO  system  for  the  purpose  of  “diagonalizing” 
the  plant  matrix.  Such  transformations  would  make  the  task  of  designing  and 
applying  a  feedback  controller  much  easier  because  the  MIMO  system  is  reduced 
to  a  set  of  decoupled  SISO  systems.  This  is  particularly  advantageous  for  systems 
with  large  numbers  of  inputs  and  outputs.  In  the  last  chapter  we  introduced 
recursive  orthogonal  transforms,  a  class  of  linear  operators  which  are  designed  to 
take  advantage  of  the  parallel  processing  capability  of  a  control  network  while 
eliminating  the  need  for  global  communication.  Here  we  begin  to  bring  these  two 
ideas  together  by  showing  how  to  find  ROTs  which  most  closely  diagonalize  a 
constant  matrix.  The  problem  of  finding  the  appropriate  ROT  variables  is  solved 
using  a  gradient  flow  which  turns  out  to  be  closely  related  to  the  so  called  double 
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bracket  equation  introduced  by  Brockett  [10]  and  discussed  in  detail  by  Helmke 
and  Moore  [22],  This  is  no  coincidence  since  our  work  is  heavily  motivated  by 
these  authors. 

Here  we  assume  that  the  ROT  parameters  are  given  and  fixed.  In  Section 
4.1  we  show  how  find  the  ROT  variables  so  that  the  resulting  ROT  most  nearly 
diagonalizes  a  given  real-valued  symmetric  matrix.  Section  4.2  describes  how  to 
find  ROT  variables  to  approximate  the  SVD  factors  of  a  given  complex-valued 
matrix. 


4.1  Diagonalization  of  Symmetric  Matrices 

The  objective  of  this  section  is  stated  as  follows:  Given  a  symmetric  n  x  n  matrix 
Ho  and  configuration  parameters  for  the  ROT 

e  =  p1o1p2o2---pLoL,  (4.i) 

find  ©i,  ©2 ,...,©l  such  that  the  matrix 

H  =  QTH0Q  (4.2) 

is,  in  some  sense,  most  closely  diagonalized.  Our  first  step  toward  solving  this 
problem  is  to  find  a  “diagonalness”  functional  < p(H ).  Then  we  search  for  the 
(©i,  @2, . . . ,  ©l)  which  minimize  by  flowing  along  the  gradient  vector  field  V(f> 
on  the  configuation  space  of  the  ©ds.  This  idea  is  motivated  by  Brockett  [10],  who 
showed  that  the  matrix  diagonalization  problem  can  be  solved  using  by  solving  a 
matrix  ODE  which  evolves  on  the  group  of  orthogonal  matrices.  As  a  result,  we 
begin  this  discussion  with  a  review  his  work.  We  then  return  to  our  problem  and 
construct  a  gradient  flow  designed  to  find  the  appropriate  ROT  variables. 
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4.1.1  Matrix  Diagonalization  via  Flows  on  SO(n ) 

In  1988,  Brockett  showed  that  a  number  of  interesting  problems  in  linear  algebra 
and  computer  science  could  be  solved  by  the  use  of  dynamical  systems  defined  on 
matrix  groups[10].  Following  [22],  we  review  his  work  on  the  diagonalization  of 
symmetric  matrices. 

Finding  the  Gradient  Vector  Field 

Let  Ho  be  a  real  valued,  symmetric,  n  x  n  matrix.  The  goal  is  to  find  0  in  SO(n ) 
such  that  H  =  @ri7o0  is  diagonal.  Note  that  H  is  symmetric  for  every  0  G  SO(n ) 
since 

Ht  =  ©Ttf0T© 

=  eTH0e 

=  H.  (4.3) 

Define 

J(0)  =  ||iV-0rilo0||2,  (4.4) 

where  N  is  a  diagonal  matrix  with  distinct  values  on  the  diagonal  and  ||-||  is  the 
Frobenius  norm,  ||A||2  =  tr  (ATA).  Hence,  J  is  the  distance  between  QTH0Q  and 
a  fixed  diagonal  matrix.  Intuitively,  minimizing  J  makes  H  =  077/(|0  look  as 
much  like  N  as  possible.  Since  N  is  diagonal,  the  0  which  minimizes  J  should 
diagonalize  H.  Simple  matrix  manipulation  reveals 

|| N  -  eTH0Q ||2  =  \\N\\2  +  ||0Ttfo0||2  -  2tr  (NH) .  (4.5) 

The  matrix  N  is  constant,  so  ||iV||2  is  constant.  The  matrix  Hq  is  constant  and  0 
is  orthogonal,  so  ||0Ti/o0||"  =  ||-^o||2  is  also  constant.  As  a  result,  the  problem  of 
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minimizing  J(0)  can  be  posed  as  the  problem  of  maximizing  the  functional 


0(0)  =  tr  ( NOTH0Q )  .  (4.6) 

Brockett’s  idea  is  to  search  for  the  maximizing  0  by  flowing  along  the  gradient 
vector  held  of  0  on  SO(n).  Recall  that  the  gradient  vector  held  of  0  on  SO(n ) 
satishes  the  following  two  properties: 

1.  V0(0)  G  TeSO(n )  for  all  0  G  SO(n). 

2.  D0e(O  =  (V0(0),O  for  all  £  G  TeSO(n). 

As  a  result  of  the  hrst  property,  V0(0)  must  be  of  the  form  0ff,  where  G  so (n). 
This  means  that  V0(0)  can  be  written  in  the  form 

V0(0)  =  0Sdv<we),  (4.7) 

where  f2vA0)  is  skew-symmetric  for  each  0  G  SO(n).  Combining  this  with  the 
second  property  gives 

L>0©(0fl)  =  (0ffv^(e),01d> 

=  — tr  (Jdv^(0)n)  .  (4.8) 

Taking  the  directional  derivative  of  0, 

T>0©(0 ft)  =  tr  ( NflTQTH0Q  +  NQTH0en)  .  (4.9) 

Using  the  skew-symmetry  of  and  the  fact  that  tr  (AB)  =  tr  (BA)  we  get 

T>0©(0Jd)  =  tr  ((NQT H0Q  —  QtH0QN)  fl) 

=  tr([N,eTH0e]n).  (4.10) 

Substituting  this  into  4.8  allows  us  to  restate  the  second  condition  as 

tr  ( [TV,  OtHo0]  n)  =  -tr  (Sdv<Ke)Sd)  .  (4.11) 


77 


So  the  choice  Dv<^0)  =  —  [ TV ,  0Tffo@]  satishes  the  second  condition.  This  is  a 
valid  choice  because  the  matrix  [TV,  C-)7  //0C-)]  is  skew-symmetric,  which  is  easily 
seen  as  follows: 

[N,H]t  =  ( NH-HN)t 

=  htnt  —  ntht 

=  HN-NH 

=  ~[ N,H ].  (4.12) 

As  a  result,  the  gradient  vector  held  is  given  by  the  expression 

V0(0)  =  -0  [TV,  eTH0e]  .  (4.13) 

The  ascent  direction  gradient  how  is  then  given  by  the  ODE 

0  =  -0  [TV,  eTH0Q]  ,  0(0)  =  0O.  (4.14) 

Convergence  Properties 

From  the  properties  of  a  gradient  how  on  a  compact  manifold  listed  in  Section  3.1.1 
it  is  clear  that  the  solution  of  Equation  4.14  exists  for  all  time  t  >  0  and  converges 
to  an  equilibrium  point  for  every  O0  G  SO(n).  Let  0*  be  such  an  equilibrium 
point  and  dehne  H  =  0(7/()0*.  Since  0*  is  nonsingular,  the  matrix  [TV,  H]  must 
be  zero.  The  (f,j)th  element  of  [TV,  if]  is 

[N,Hlu  =  (NH)lt  -  (HN)V 

n  n 

=  NjkHkj  —  HjkNkj.  (4-15) 

k= 1  k= 1 

The  matrix  TV  is  diagonal  so  [TV]^  =  0  for  i  ^  k.  Let  rq  denote  the  Tth  diagonal 
element  of  TV  and  let  denote  the  (i,j) th  element  of  H.  Then 

[TV,  H\i-  =  riihij  -  hijrij 
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hij  ( Hi  —  VLj)  . 


(4.16) 


The  diagonal  elements  of  N  are  distinct,  so  for  i  ^  j,  [. N ,  H]^  =  0  only  if  =  0. 
Therefore  the  matrix  H  must  be  diagonal. 

4.1.2  Matrix  Diagonalizing  Flows  for  ROTs 

Here  we  extend  the  results  of  Brockett  [10]  to  find  the  recursive  orthogonal  trans¬ 
form  0  which  most  closely  diagonalizes  a  real-valued  symmetric  matrix. 

Preliminaries  and  Definitions 

Let  H0  be  a  real  valued,  symmetric,  n  x  n  matrix.  This  is  the  matrix  we  seek  to 
diagonalize. 

Let 


d\ 

0 

...  0 

— 

0 

0k2 

...  0 

,  (4.17) 

0 

0 

for  k  —  1, 2, . . . ,  L,  some  L  >  0,  where  9j  G  SO(rikj), 

mk 

^^rikj  =  n  Vfc  =  1, 2, . . . ,  L,  (4-18) 

5= 1 

and  the  Os  represent  appropriately  dimensioned  blocks  of  zeros.  In  simple  language, 
each  0 is  a  block  diagonal  matrix  where  the  blocks  on  the  diagonal  are  orthogonal 
matrices.  Let  Mk  =  SO(nk i)  x  SO{nk 2)  x  •  •  •  x  SO{nkmk)-  Then  0^  belongs  to 
the  Lie  subgroup  Mk  C  SO(n). 

Recall  that  the  tangent  space  of  SO(£)  at  the  identity  is  the  Lie  algebra 

so (i)  =  {fl£  R£x£|flT  =  -fi}.  (4.19) 
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The  Lie  algebra  of  Mk  is  then 


TeMk  =  so(nki)  ®  so (nk2)  ®  •  •  •  ©  so (nfemJ.  (4.20) 

A  vector  in  TeMk  can  then  be  written 

0  •  •  •  0 

0  w!  -  0 

■  (4.21) 

_  0  ...  0  _ 

where  tVj  G  so(nkj )  for  j  =  1,2, .. . ,  m*,.  Hence  the  tangent  space  to  Mk  at  the 
point  0fc  is 

TekMk  =  {QkQk  ^  TeMk}  .  (4.22) 

Dehne  T  G  M  =  Mi  x  M2  x  •  •  •  x  Ml  to  be  the  ordered  L-tuple  of  ©^s, 

T  =  (01,02,...,0l).  (4.23) 

Likewise,  let  X  denote  a  vector  in  Tq,M, 

x  =  (e1n1,e2n2,...,eLnL).  (4.24) 


For  k  =  1,  2,  3, . . . ,  L,  Pk  is  a  fixed  permutation  of  the  n  x  n  identity  matrix. 
Each  Pk  is  an  element  of  the  group  0  (n) . 

Finally,  for  k  =  1,2 , . . . ,  L,  we  dehne  two  sequences  of  recursive  orthogonal 
transforms: 

k 

Qk  =  (4.25) 

£= 1 
n 

=  n p^-  (4-26) 

e=k+i 
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Here  the  product  symbol  denotes  multiplication  with  indices  ascending  from  left 
to  right.  Also,  \]^k=a  =  I  for  b  <  a.  We  call  ©*,  the  fcth  lower  subtransform  and 
(-) /.  the  kth  upper  subtransform.  The  ROT  0  =  0  L  is  referred  to  as  the  complete 
transform. 

Constructing  the  Gradient  Flow 

The  objective  of  this  section  is  to  approximately  diagonalize  the  symmetric  matrix 
Hq.  This  is  done  by  finding  T  =  (0i,  02, . . . ,  ©l)  such  that  the  matrix 

H(V)  =  QtlP[  ■  ■  •  0^P2r0fP1rJHoPi0iP202  •  •  •  PlQl 

=  eTH0e  (4.27) 

is  “more  diagonal”  than  it  is  for  any  other  choice  of  T. 

Towards  this  goal,  let  A"  be  a  fixed  diagonal  n  x  n  matrix.  We  can  now  define 
a  cost  function  using  distance  between  H  and  N  given  by  the  Frobenius  norm. 
Simple  matrix  manipulation  reveals 

|| N  -  H(V) ||2  =  ||iV||2  +  ||tf  (^)f  -  2tr  (NH(V)) .  (4.28) 

Since  N  is  hxed,  ||  W||2  is  constant.  Since  Ho  is  hxed  and  the  matrix  0  is  orthogonal, 
~  ~  2 

||i7||  =  0Ti7o©  is  also  constant.  The  Frobenius  distance  between  H  and  N  is 
minimized  when 

0(T)  =  tr(iVW(T))  (4.29) 

is  maximized.  The  function  (f>  :  M  — »  H  can  be  thought  of  as  a  “diagonalness” 
functional. 

The  idea  here  is  to  search  for  the  T  which  maximizes  (j)  by  flowing  along  the 
gradient  vector  held  of  (j)  on  the  manifold  M.  Recall  that  V0(T)  on  M  is  defined 
by  the  following  properties: 
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1.  V0O)  e  TyM  VT  e  M. 


2.  B0*(x)  =  (V0(^), x)  vx  e  t*m. 

The  first  property  states  that  V0(T)  in  contained  in  the  tangent  space  of  M  at 
the  point  \F  This  means  that  V0(T)  must  be  of  the  form 

V0(T)  =  (©iQ^, . . . ,  ©l^I0)  ,  (4.30) 

where  G  TeMk ,  A:  =  1,2,  For  convenience  of  notation,  we  have  sur- 

pressed  the  fact  that  depends  on  T. 

From  the  second  property,  we  must  have 


£>^((©1^!,  . . . ,  ©ifd^))  —  (V0(©i,  . . . ,  ©l),  (©i^i,  . . . ,  ©l^l)) 

=  ((0!^,  •  •  •  ,  ©L^),  (©1^1,  •  •  •  ,  ©L^l)} 

=  J^tr  ((Sd^T5^)  .  (4.31) 

k= 1 


Finding  an  expression  for  Dcp^^X): 

dmx)  =  tr n  p&i\\ 


^=1 


i=i 


Ki=l 


'k- 1 


e=k+ 1 


+iv^n^j  ^o^n^0<j^©A  ( n 

/  L  _  T 

tr  ^jv  (eknkek^  b0©  +  NeTH0eknkek 

\k= 1 

tr  ^7V©^©^o©  +  AV©T^o©A©a 


(4.32) 


V  /c=l 


Using  tr  (AB)  =  tr  (BA)  and  the  skew-symmetry  of  Qk  we  get 


Dfa(X)  =  tr  (  J2  (®kNOTH0Ok  -  ©^BO07V©^  Q 

,fc=i 


(4.33) 
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Using  the  fact  that  0  =  0^0^  for  any  k  —  1,  2, . . . ,  L,  we  get 


DMX)  =  tr  ^(0fciV0j0^o0fc-0^o0fc0fciV0r)a 

\k= 1 

\k= 1 

For  k  —  1, 2, . . . ,  L,  define 

Hk  =  eTkH0ek 
Nk  =  ekNQl. 


(4.34) 


(4.35) 

(4.36) 


The  matrices  Hk  and  Nk  are  symmetric  for  each  k.  This  means  that  the  bracket 
[Nk,Hk]  G  TeO{n).  In  order  to  get  a  valid  gradient  flow,  each  Qk  in  Equation  4.34 
must  be  multiplied  on  the  left  by  (f!^)T,  where  G  Te(Mk).  To  accomplish 
this,  we  introduce  a  set  of  natural  projection  operators,  IIfc  :  TeO{n )  — *  Te(Mk). 
IU  projects  from  the  set  of  skew-symmetric  matrices  to  the  set  of  block  diagonal 
skew-symmetric  matrices  by  setting  all  off-diagonal  blocks  to  zero. 

Fact  4.1.1  Let  A  G  o  (n)  and  let  Lli  G  TeMi  as  defined  in  Equation  f.21.  Then 

tr  ((Iljkl)Tflj)  =  tr(ATClf)  .  (4.37) 


Proof: 

Looking  at  Equation  4.21  we  see  that  the  (j,  fc)th  block  of  fh  is 

{fli\jk  =  u Jij53k,  (4.38) 

where  8k  is  the  Kronecker  symbol.  Partition  A  into  blocks  a jk  so  that  the  block 
structure  of  A  is  compatible  with  that  of  fl*.  Then 

rrii  rrii 

tr  (. ATQl )  =  EE  tr(a jk[^i\jk)  (4.39) 

j= 1  k= 1 
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(4.40) 


rrii  mi 

=  EEtr  (a 

j= 1  fc=l 

mi 

^  ]  tr 

fc=i 

=  tr((n^)Ta). 


(4.41) 

(4.42) 


In  light  of  this  fact, 

L 

D<h(X )  =  £tr((-nfe  [iVfc,hffc])Ta)  •  (4.43) 

fc=i 

If  we  let  =  —11*.  [iVfc,#fc],  then  V0(\l/)  meets  the  conditions  for  a  gradient 
vector  held  outlined  above.  This  yields 

V0(T)  =  (-©in,  [TV1,#1]  ,-02n2  [TV2,#2]  ,...,-eL nL  [tvl,#l])  .  (4.44) 


The  corresponding  ascent  direction  gradient  how  is  then  given  by  a  system  of 
coupled  matrix  ODEs.  This  system  is  written  as 


QkNQl,QlHnQ, 


©fc(0)  —  ©fco, 


(4.45) 


for  k  =  1, 2, . . . ,  L. 

At  this  point  it  is  natural  to  bring  up  the  question  of  the  convergence  properties 
of  this  gradient  how.  From  the  properties  of  gradient  hows  on  compact  manifolds 
listed  in  Section  3.1.1  we  know  that  the  solution  to  Equation  4.45  exists  for  all  time 
and  converges  to  an  equilibrium  point  for  every  set  of  initial  conditions.  Except 
in  a  few  special  cases  (see  Appendix  A)  it  is  difficult  to  say  more  than  this.  The 
numerical  results  we  have  obtained  are  promising,  however,  as  is  demonstrated  by 
the  following  example. 
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Example 


Here  we  show  the  results  of  applying  this  technique  to  approximately  diagonalize  a 
16  x  16  symmetric  matrix  with  random  entries.  Here  we  used  the  ROT  structure  for 
a  ID  array  of  sensor  actuator  pairs  presented  in  Section  3.2.1.  The  original  plant 
matrix  is  shown  in  Figure  4.1.  The  approximately  diagonalized  plant  matrix  for 
ROT  levels  3,7,11,  and  15  are  shown  in  Figures  4.2,  4.3,  4.4,  and  4.5,  respectively. 
Figure  4.6  show  a  plot  of  approximation  error  as  a  function  of  the  number  of  levels 
used  in  the  approximating  ROT.  Here,  the  approximation  error  is  defined  to  be 


E(Q) 


QTGQ  -  0TG0 

INI 


(4.46) 


where  Q  is  a  matrix  whose  columns  are  the  unit  eigenvectors  of  G. 

The  dimension  of  SO[n )  is  n(n— 1)/2.  The  number  of  degrees  of  freedom  in  the 
ROTs  being  considered  is  Ln/ 2,  where  L  is  the  level  of  the  transform.  Intuitively, 
when  L  =  n—  1  the  ROT  “should”  have  enough  degrees  of  freedom  to  represent  any 
0  in  SO(n).  These  results  seem  to  support  this  notion  since  the  approximation 
error  gets  very  close  to  zero  for  L  =  15. 


4.2  Singular  Value  Decomposition 

In  Section  4.1  we  addressed  the  problem  of  finding  a  recursive  orthogonal  transform 
to  diagonalize  a  real-valued  symmetric  matrix.  In  this  section,  we  extend  these 
results  to  find  approximations  of  the  SVD  factors  which  diagonalize  a  general 
complex-valued  matrix.  In  this  case,  given  an  nxm  complex  matrix  //(l  and  ROT 
configurations  for  U  G  U(m), 

Lu 

U  =  l[PiUi:  (4.47) 

i=  1 
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and  V  G  U(n), 


Ly 

V  =  Y[QiVi,  (4.48) 

i= 1 

we  will  search  for  the  block  diagonal  Uts  and  VjS  so  that  the  matrix  H  G  (Drixri, 

H  =  VhH0U  (4.49) 

is  most  closely  diagonalized.  Recall  that,  according  to  our  definition,  a  matrix  A 
is  diagonal  if  [A\i3  =  0  for  all  i  ^  j.  This  definition  applies  to  non-square  matrices 
as  well  as  square  ones. 

As  in  the  case  of  diagonalizing  symmetric  matrices,  our  approach  is  to  search  for 
the  best  UiS  and  Rs  by  flowing  along  the  gradient  vector  held  of  a  “diagonalness” 
functional.  Here  we  are  motivated  by  Helmke  and  Moore  [21],  who  extended 
Brockett’s  work  on  symmetric  matrices  to  perform  SVD  via  a  gradient  how  on 
U{m )  x  U{n).  We  begin  this  discussion  with  a  review  of  their  work.  We  then 
return  to  the  problem  of  SVD  using  recursive  orthogonal  transforms. 

4.2.1  SVD  via  Flows  on  U(m)  x  U(n) 

In  1992,  Helmke  and  Moore  [21]  addressed  the  problem  of  finding  the  SVD  fac¬ 
tors  of  a  complex-valued  matrix  by  using  gradient  hows  on  the  unitary  group. 
This  work  is  an  extension  to  the  work  or  Brockett  [10],  which  we  reviewed  in 
Section  4.1.1.  Here  we  present  a  brief  summary  of  the  work  in  [21]. 

Finding  the  Gradient  Flows 

Given  H0,  N  G  (Dnxm,  consider  the  task  of  finding  U  G  U(m ),  V  G  U(n)  to  minimize 
the  squared  Frobenius  distance 

|| N  -  VhH0U\\2  =  ||iV||2  +  ||Ho||2  -  2retr  ( NhVhH0U )  .  (4.50) 
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N  and  Ila  are  constant,  so  minimizing  the  above  norm  is  equivalent  to  maximizing 
the  functional 

cj)(U,  V )  =  2retr  ( NhVhH0U )  .  (4.51) 

If  the  matrix  N  is  diagonal,  then  <f>  can  be  thought  of  as  a  “diagonalness”  function 
for  the  matrix  VH HQU . 

As  in  the  case  of  Brockett’s  work,  the  idea  is  to  search  for  the  U  and  V  which 
maximize  (f>  by  flowing  along  the  gradient  vector  held  of  0  on  the  smooth  Lie  group 
U(m)  x  U[n). 

The  Lie  algebra  of  U{m )  x  U(n)  is  appropriately  written  as  the  direct  sum 
u(m)  ©  u(n).  Accordingly,  a  tangent  vector  at  the  point  (U,  V )  is  a  vector  of  the 
form  (LT2,  VA),  where  G  u(m)  and  A  G  u(n). 

Recall  the  properites  of  the  gradient  how: 

1.  V(I>(U,  V )  G  T(uy)  ( U(m )  x  U(n))  for  all  {U,  V )  G  U(m)  x  U(n). 

2.  D<j){uy)(im,  VA)  =  V ),  {im,  VA))  for  all 

(un,v. A)  G  T(uy)  (U(m)  x  U (n)). 

The  hrst  property  says  that  V0(C/,  V)  is  tangent  to  U(m)  x  C/(n)  at  every  point, 
which  means  that  V(j)  can  be  written 

V<j)(U,  V )  =  (£/$dV0,  VAV(t>)  ,  (4.52) 

where  (fiV9i,  Av'?i)  G  u (m)  ®u(n).  For  convenience  of  notation,  we  have  surpressed 
the  fact  that  and  Av ^  depend  on  ([/,  V).  Writing  V0  in  this  way,  the  second 
property  becomes 

D<f>{UtV)(Un,V  A)  =  ((f/flv^,RAv^),(175d,l/A)> 

=  — retr  —  retr  (Av?iA)  .  (4.53) 
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Finding  an  expression  for  D4>{u,v){UVL,  V K)\ 


D(t){u>v)(im,V A)  =  2retr  (NhAhVhH0U  +  NHVHH0im) 

=  2retr  (. NhVhH0UVL  -  VhH0UNhA )  (4.54) 

According  to  Fact  3.1.3,  NHVH H0U  and  VH H0UNH  can  be  replaced  with  their 
skew-hermetian  components  yielding 

D<j>{uy)(im,V. A)  =  2retr  Q  (NhVhH0U  -  UhH*VN)  ft 

(1 VhHqUNh  -  NUhH^V)  A^ 

=  retr  ({TV,  ft  +  {TVH,  A)  ,  (4.55) 

where 

{•,•}  :(Dpxs  x(DpX9  -G  u(q) 

(A,B)  ^  AhB~BhA  (4.56) 

is  the  extended  Lie  bracket.  Noting  that  {TV,  VH H0U }  G  u (to)  and 

{TVH,  e  u(n)  and  looking  back  at  Equation  4.53,  we  see  that  both 

properties  of  the  gradient  vector  field  are  satisfied  for  the  choices 


ftv7  =  _  {jv,  VhH0U} 

Av^  =  -  {Nh,  UhH^V}  . 

(4.57) 

(4.58) 

This  results 

in  the  gradient  vector  held 

V<f>(U,  V )  =  ( -U  {TV,  VhH0U}  ,  -V  {Nh,  UhHqV})  . 

(4.59) 

The  resulting  gradient  flows  are  then  given  by  the  coupled  pair  of  matrix  ODEs 

evolving  on  U{m )  x  U{n).  They  are 

U  =  -  {TV,  VhH0U}  ,  7/(0)  =  U0E  U (to) 

E  =  -  {NH,  UhHqV }  ,  y(0)  =  E0  e  U(n). 

(4.60) 

(4.6!) 
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4.2.2  SVD  Flows  for  ROTs 


Here  we  extend  the  work  of  Helmke  and  Moore  [21]  to  find  the  recursive  orthogonal 
transforms  U  and  V  which  most  closely  approximate  the  SVD  factors  of  a  given 
matrix  Hq. 


Preliminaries  and  Definitions 

Let  H0  G  (Dnxm  be  the  matrix  to  be  decomposed.  Let  N  be  a  fixed  matrix  in  07 
Let 

r 

u\  0  ■  ■  ■  0 

•  0 


Ui  = 


0  u\ 


0  ■  ■  ■  0  ul 


for  i  —  1, 2, . . . ,  Li /,  some  Lv  >  0,  where  ui-  G  U(rriij), 


(4.62) 


ki 


y:  rriij  =  m  Vi  =  1, 2, . . . ,  Ln, 

3= 1 

and  the  Os  represent  appropriately  dimensioned  blocks  of  zeros.  Let 


(4.63) 


Mf  =  U(mu )  x  U(mi2)  x  •  •  •  x  U(miki). 


(4.64) 


Then  Ui  belongs  to  the  smooth,  connected  Lie  subgroup  C  U(m).  The  Lie 
algebra  of  Mf  is  then 

TeMV  =  u (ma)  ©  u (mi2)  ©  •  •  •  ©  u (miki).  (4.65) 


We  choose  to  write  and  element  of  TeMf  as  a  block  diagonal  skew-hermetian  with 
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the  same  block  structure  as  Ui,  i.e.  G  TeM^r  will  be  written 


a 


u[  0  •  •  •  0 
0  u>2  ■  ■  ■  0 


o  o  u4 


(4.66) 


where  a;*  G  u{rriij)  for  j  =  1,2 fa.  The  tangent  space  of  M J7  at  the  point  Ui  is 
then 


TUtM-J  =  {U,nt  1 0,;  G  TeM?  }  . 
We  make  a  similar  set  of  definitions  for  Vi.  Let 


v\  0  ■  ■  ■  0 
0  v\  ■  ■  ■  0 


0  ■  ■  ■  0  v\. 


for  i  —  1, 2, . . . ,  Ly ,  some  Ly  >  0,  where  u*  G  U(riij ), 


h 

=  n  Vi  =  1, 2, ... ,  Lv. 
j= i 


(4.67) 


(4.68) 


(4.69) 


Let 


MY  =  U(na )  x  £7(71*2)  x  •  •  •  x  U(nUi).  (4.70) 

Then  V)  belongs  to  the  smooth,  connected  Lie  subgroup  MY  C  U(n).  The  Lie 
algebra  of  MY  is 


TeM Y  =  u(nn )  ®  u(ni2)  ©  •  •  •  ®  u(n^).  (4-71) 


As  before,  we  choose  to  write  and  element  of  TeMY  as  a  block  diagonal  skew- 
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hermetian  with  the  same  block  structure  as  Vt .  i.e.  Aj  G  TeM. Y  will  be  written 


X[  0  •  •  •  0 
0  A^  • • •  0 


0  0  Aj. 


(4.72) 


where  A*  G  u (n^)  for  j  =  1,2,...  ,  £j.  The  tangent  space  of  at  the  point  Vj  is 
then 

TViMY  =  {ViK  \K  G  TeMY  }  .  (4.73) 


Define  the  smooth  Lie  groups 


Mu 

=  Mf  x  M%  x  • 

•  •  x  mV 

Lu 

Mv 

=  M\  x  Ml  x  • 

■■xMY 

LV 

M 

=  Mu  x  My. 

(4.74) 


Let  f  G  M  be  an  (. Ljj  +  Ly)-tuple  of  Uts  and  VjS, 


*  =  (U1,U2,...,ULu,V1,V2,...,  VLv) .  (4.75) 

Then  X  G  T^M  is  an  (Lrj  +  Ly  )-tuple 


x  =  {u^,  u2n2, . . . ,  uLunLu,v1  au  v2a2 , . . . ,  vLvaLv)  ,  (4.76) 

where  0*  G  TeM^ ,  7  =  1,2,...,  Ljj ,  and  Aj  G  TeM^ ,  i  =  1,2,...,  Ly. 

The  matrices  Pi,  i  =  1,2,...,  Ljj,  are  fixed  permutations  of  the  m  x  m  iden¬ 
tity  matrix.  Likewise,  Qj,  i  =  l,2,...,Ly,  be  hxed  permutations  of  the  nxn 
identity  matrix.  Both  sets  of  matrices  are  given  as  part  of  the  ROT  configuration 
parameters. 
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For  i  =  1,2, ,  Lu,  we  define  the  sequences  of  rn  dimensional  ROTs 

i 

U,  =  n  PjU,  (4.77) 

3  = 1 
Lu 

Ui  =  n  piui-  <4-78) 

j=i+ 1 

Likewise,  for  i  =  1,2, ,  Ly  we  define  the  n  dimensional  ROTs 

i 

^  =  n«(»i  <4-7s)) 

i=i 

Lv 

K  -  n  0^..  (4,80) 

J=*+l 

We  will  call  these  ROTs  subtransforms;  Ui  and  Id  are  lower  subtransforms,  while 
Ui  and  Vi  are  upper  subtransforms.  The  ROTs  U  =  Ulv  and  V  =  V}  v  will  be 
called  complete  transforms. 


Constructing  the  Gradient  Flow 

The  objective  of  this  section  is  to  find  recursive  orthogonal  transforms  U  and  V  of 
the  configuration  outlined  above  that  most  closely  approximate  the  SVD  factors  of 


the  matrix  H0.  To  do  this,  we  will  search  M  for  the  T  =  (U\, . . . ,  Ulv,  Vi, ,  Vlv) 
which  most  closely  diagonalizes  the  matrix 

H(V)  =  VhH0U.  (4.81) 

Consider  the  Frobenius  distance  between  the  fixed  matrix  N  and  II (T), 

|| N  -  H(V) ||2  =  ||iV||2  +  ||tf0||2  -  2retr  (NhH(V))  .  (4.82) 

Clearly,  this  quantity  is  minimized  when  the  functional 

=  2retr  (NhH($>))  (4.83) 
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is  maximized.  If  we  let  the  matrix  IV  be  a  diagonal  matrix  with  distinct  values  on 
the  diagonal,  then  0  can  be  thought  of  as  a  “diagonalness”  function. 

As  in  the  case  of  diagonalizing  symmetric  matrices,  the  approach  we  take  is  to 
search  the  Lie  group  M  for  the  maximizing  T  by  flowing  along  the  gradient  vector 
held  of  (j).  The  gradient  vector  held  V0  is  dehned  by  the  following  properties: 


1.  V0(T)  G  T^M  VT  G  M. 

2.  Dfo(X)  =  (V<f>{V),X)  VX  G  T*M. 


To  satisfy  the  hrst  property,  V0(\&)  must  be  of  the  form 


V0(T)  =  (V,  si**,.. 


(4.84) 


where  G  TeM^ ,  i  =  1,2  and  G  TeM^ ,  i  =  1,2,  ...,Ly.  For 

convenience  of  notation,  we  have  surpressed  the  fact  that  these  matrices  depend 
on  T.  Writing  V0  in  this  form,  the  second  property  states 


Lu  Ly 

=  -  retl  ~  retr 


X 


i= 1  i= 1 

Finding  and  expression  for  D(f><s,(X),  we  get 

Ljj  Ly 


l  _  _  v  _  /  H 

DMX)  =  2Tetrlj2NHVHHoUiniiji  +  J2NH{ViKV?)  H0U 


v  i=l 
Lu 


i= 1 
Lv 


2retr  ^  UiNHVH H 0[/A  -  ^  V,H H0UNHVtH ^ 


K  i= 1 
(  Lu 


i— 1 


\  i—i 


2retr  ^  -  ( UiNHVHH0Ux \  -  U^H^VNU?) 


1=1 


(4.85) 


-  XI  F  (v?H0UNHV?  -  VtNUHH«V<)  A 
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=  J>etr  ({ViNU^VfHoUi}^ 

i= 1 

+  ^retr  ({ U.A^Vf ,  V? A,)  (4.86) 

i= 1 

In  the  last  Equation,  the  matrices  multiplying  0,:  and  A,:  are  in  u (to)  and  u(n), 

respectively.  Recall  that  in  order  for  the  gradient  vector  to  be  tangent  to  M, 

each  0,;  must  be  multiplied  on  the  left  by  Jiff  where  G  Te(MP)  and  each 

A i  must  be  multiplied  on  the  left  by  Aff  where  Af^  G  Te(Mf).  To  accomplish 

this,  we  introduce  a  set  of  natural  projection  operators,  Ilf  :  TeU(m )  — *  Te(Mf ), 

i  =  1, 2, . . . ,  Lu  and  Ilf  :  TeU(n )  — »  Te(M ]  ),  i  =  1, 2, ... ,  Ly.  These  operators 

project  from  the  set  of  skew-hermetian  matrices  to  the  set  of  block  diagonal  skew- 

hermetian  matrices  by  setting  all  off  diagonal  blocks  to  zero.  Using  these  operators 

along  with  Fact  4.1.1,  we  get 
Lu 

DM* 0  =  X>tr(nf  {ViNU^VfHoUi}^ 

1=1 

+  X>tr  (nr  {UiNHVf,  Up Hq  Vi }  A,,)  .(4.87) 
1=1 

Looking  back  to  Equation  4.53,  we  see  that  the  choices 

^  =  -nf  {ViNUPiVfHoUi},  i  =  l,2,...Lu, 

A ?+  =  -n Y  {uNhVP,UPhPV^  i  =  1,2,...Lv.  (4.88) 

satisfy  both  gradient  flow  properties.  As  a  result,  the  gradient  flow  equations  are 
Ui  =  -UiUY  [ViNUP,  VPH0U}  UP 0)  =  Ui0  G  Mf ,  (4.89) 

for  i  =  1, 2, . . . ,  Lv,  and 

Vi  =  -VJlY  {uNhvP,  UpHPVi }  ,  vpo)  =  Vi0  G  mV,  (4.90) 

for  i  —  1, 2, . . . ,  Ly. 
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Plant  Matrix  (1 D  Membrane) 
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Figure  4.7:  Plant  matrix  for  flexible  membrane. 

Examples 

Here  we  apply  this  method  to  the  examples  introduced  in  Section  2.4.  The  plant 
matrix  for  the  simplified  membrane  example  is  depicted  in  Figure  4.7.  Approxi¬ 
mately  diagonalized  plant  matrices  using  10,  20,  and  31  level  ROTs  are  shown  in 
Figures  4.8,  4.9,  and  4.10,  respectively.  The  plant  matrix  for  the  flexible  beam 
at  its  6th  resonance  is  plotted  in  Figure  4.11.  Approximately  diagonalized  plant 
matrices  using  10,  20,  and  31  level  ROTs  are  shown  in  Figures  4.12,  4.13,  and  4.14, 
respectively.  These  results  will  be  discussed  more  quantitatively  in  Chapter  6. 
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Chapter  5 


Spatially  Invariant  Systems 


In  many  applications  of  sensor /actuator  arrays,  a  large  number  of  identical  sen¬ 
sor/actuator  pairs  are  placed  at  regular  intervals  in  some  homogeneous  medium. 
Such  systems  often  exhibit  a  spatial  invariance  property;  the  coupling  between 
any  two  nodes  depends  only  on  the  relative  positions  of  the  nodes  and  not  on  their 
absolute  positions  in  the  array.  This  property,  which  will  be  defined  more  rigor¬ 
ously  below,  can  be  exploited  to  aid  in  the  design  and  implementation  of  feedback 
controllers.  Brockett  and  Willems  [11]  [12]  used  the  spatial  invariance  property 
found  in  discretized  partial  differential  equations  to  develop  optimal  feedback  laws. 
Melzer  and  Kuo  [27]  formulated  a  Riccatti  equation  for  spatially  invariant  systems 
and  then  applied  their  results  to  the  problem  of  controlling  an  infinite  string  of 
vehicles.  El-Sayed  and  Krishnaprasad  [16]  employed  a  similar  method  to  design 
optimal  feedback  laws  to  control  the  depth  of  a  seismic  cable  containing  an  array 
of  hydrophones.  More  recently,  Bamieh  [1]  and  Bamieh,  Paganini,  and  Dahleh 
[2]  have  proposed  the  use  of  this  same  idea  to  develop  controllers  for  large  scale 
sensor /actuator  arrays. 

Here  we  propose  to  exploit  the  spatial  invariance  property  for  the  purpose  of 
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plant  matrix  diagonalization.  We  begin  by  defining  a  spatially  invariant  system 
and  discussing  some  of  the  features  that  such  systems  exhibit.  We  then  show  that 
the  discrete  Fourier  transform  can  be  used  to  diagonalize  the  plant  matrix  of  any 
spatially  invariant  system.  Using  this  information,  we  then  address  the  problem 
of  diagonalizing  such  systems  using  ROTs. 

5.1  Infinite  Sensor /Actuator  Arrays 

Here  we  define  the  spatial  invariance  property  for  an  infinite  sensor/actuator  array 
and  we  discuss  how  this  property  can  be  exploited  to  aid  in  controller  design  and 
implementation  for  such  systems.  Of  course,  this  discussion  is  purely  hypothetical 
since  infinite  arrays  do  not  exist.  In  the  next  section  we  will  extend  the  ideas 
discussed  to  arrays  of  finite  length. 

Consider  a  control  system  consisting  of  some  physical  plant  equipped  with  a 
doubly  infinite  sensor/actuator  array.  Assume  that  each  node  on  the  array  contains 
exactly  one  sensor  and  one  actuator.  Let  yt  be  the  output  measured  at  the  zth 
sensor  and  let  Ui  be  the  input  applied  to  the  fill  actuator.  This  system  is  spatially 
invariant  if,  for  each  i,  the  output  y*  can  be  written  as 

Vi=  ■»  C5-1) 

i—jeM 

where  gr-}  is  a  scalar,  finite  dimensional  linear  operator  representing  the  dynamic 
coupling  between  yi  and  and  u3.  The  set  M  C  %  defines  the  extent  of  the  spatial 
coupling  between  y  and  u.  The  coupling  function  g^j  is  zero  whenever  i  —  j 
is  outside  of  M.  For  example,  nearest  neighbor  coupling  can  be  represented  by 
letting  M  =  {—1,0,1}.  This  system  is  invariant  under  discrete  translations.  In 
other  words,  the  relationship  between  y \  and  Uj  is  identical  to  the  relationship 
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between  yi+k  and  Uj+k  for  all  feel 

In  the  time  domain,  g^3  represents  the  convolution  (in  time)  with  a  weighting 
function  wt_3{t).  Equation  5.1  is  then  be  written  in  the  time  domain  as 

Vi(t)=  ^2  /  -  a)uj(a)da.  (5.2) 

i-j£M''0 

In  the  frequency  domain,  gt^3  represents  multiplication  by  a  hnite  dimensional 
transfer  function,  gi_3{s).  In  this  case,  Equation  5.1  becomes 

Vi(s )  =  3i-j(s)u3(s),  (5.3) 

i-j£M 

where  yi(s)  and  Uj(s)  are  the  Laplace  transforms  of  the  time  varying  quantities 
yi(t)  and  Uj(t),  respectively.  For  convenience,  the  remainder  of  this  discussion  will 
use  the  s  domain  representation  given  by  Equation  5.3. 

It  is  worthwhile  at  this  point  to  take  a  close  look  at  the  variables  in  Equation  5.3. 
First  consider  the  familiar  Laplace  variable  s.  The  s  domain  is  associated  with  the 
temporal  frequency  domain  by  making  the  assignment  s  =  jujt.  Traditionally, 
the  “temporal  frequency  domain”  is  simply  called  the  “frequency  domain” ,  but  in 
this  case  it  is  important  to  specify  that  we  mean  frequency  in  time  as  opposed  to 
frequency  in  space.  The  reason  for  this  distinction  will  soon  become  clear.  The 
second  variable  in  Equation  5.3  consists  of  the  subscript  i,  i  —  j,  or  j  in  the  case 
of  y,  g ,  or  u,  respectively.  This  variable  represents  position  in  space.  Hence,  the 
subscript  describes  the  spatial  behavior  of  the  system  while  the  Laplace  variable 
describes  the  temporal  behavior. 

The  spatial  invariance  exhibited  by  Equation  5.3  is  analogous  to  the  well  known 
time  invariance  property  exhibited  by  discrete  linear  time  invariant  systems.  Here, 

g(s)  =  {•••>0-i(s),0o(s),0i(s),...} 


105 


can  be  thought  of  as  a  spatial  impulse  response.  The  elements  of  this  impulse 
response  are  transfer  functions.  The  output  y(s)  =  {. . . ,  t/_i(s),  yo(s),  yi(s), . . .} 
can  be  written  as  y(s)  =  g(s)  *  u(s),  where  the  *  operator  represents  convolution 
in  the  spatial  domain  and  u(s)  =  {. . . ,  U-i(s),  Uo(s),  U\(s), . . .}.  It  follows  then 
that  we  can  write 

Y(s,  z)  =  G(s,  z)U(s,  z),  (5.4) 


where  Y(s,z),  G(s,z ),  and  U(s,z )  are  the  2-sided  spatial  z  transforms  of  y(s), 
g(s),  and  u(s),  respectively.  The  2-sided  spatial  z  transform  is  dehned  to  be 

OO 

X(s,z)=  Xi(s)f  (5.5) 

i— — oo 

where  i  indexes  the  spatial  domain.  Here,  z  is  the  spatial  shift  operator,  i.e. 


2  :  Xi  xi+1. 

Assume  that  y(s),  g(s),  and  u(s)  satisfy  conditions  such  that  their  z  transforms 
exist  on  the  unit  circle.  Then  the  z  domain  can  be  though  of  as  the  spatial 
frequency  domain  by  assigning  z  =  e>Wa.  Hence,  Equation  5.4  is  a  frequency 
domain  representation  of  the  system  in  both  time  and  space. 

This  representation  provides  some  advantages  in  the  development  and  imple¬ 
mentation  of  feedback  control  laws  for  the  system.  Suppose  we  wish  to  develop 
a  spatially  invariant  feedback  law  U(s,z )  =  K(s,z)Y(s,z).  Under  this  feedback, 
the  closed  loop  transfer  function  becomes 


Gd(s,z)  = 


G(s,  z) 


(5.6) 


1  +  K(s,  z)G(s,  z ) 

It  is  possible  to  evaluate  the  stability  of  the  Gci(s,z)  using  extensions  of  classical 
frequency  domain  methods  such  as  the  Nyquist  and  circle  criterion.  Under  some 
conditions,  spectral  factorization  techniques  can  be  used  to  develop  optimal  spa¬ 
tially  invariant  control  laws  [12].  Once  a  controller  is  decided  upon,  implementation 
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on  a  control  network  is  straight  forward;  the  control  law  U(s,z )  =  K(s,z)Y(s,z) 
is  applied  by  each  node  on  the  network. 

5.2  Periodic  Arrays  and  Circulant  Matrices 

In  the  previous  section  we  discussed  the  spatial  invariance  property  for  an  infinite 
sensor/actuator  array.  As  is  done  frequently  in  the  time  invariant  case,  these  ideas 
can  be  applied  to  a  finite  array  by  assuming  that  the  array  is  periodic.  Consider 
a  sensor/actuator  array  with  n  nodes  where  each  node  has  exactly  one  sensor 
and  one  actuator.  The  periodic  assumption  means  that  yi+fcn(s )  =  yi(s)  and 
Ui+kn(s)  =  Ui(s)  for  all  i  =  0,l,...,n  —  1  and  for  all  tel.  We  call  the  system 
spatially  invariant  if 

n—1 

Vi(s)  'y  ^  9((j— j)  mod  n) (s)Uj (^)  (5-7) 

3= 0 

for  i  =  0,1....,  n  —  1.  Physically,  this  means  that  the  array  must  wrap  around  so 
that  the  ( n  —  l)th  node  is  adjacent  to  the  zeroth  node.  An  example  of  this  would 
be  a  circular  array. 

Let  y(s)  =  [y0(s),yi(s), . . .  ,yn~i(s)}T  and  u(s)  =  [u0(s) ,  u^s) : ,  un_1(s)]T . 
Then  Equation  5.7  becomes 

y(s)  =  P(s)u(s )  (5.8) 

where 

,9o (s)  gn-i(s )  gn-2(s)  •••  gi(s) 

9i(s)  9o(s)  gn-i(s)  •••  g2(s) 

P(s)  =  g2(s)  (s)  5-0(5)  •••  .9.3(5)  •  (5.9) 

gn-i(s)  gn-2(s )  gn-3(s)  •••  5-0 (s) 
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In  this  equation,  gt(s)  is  a  finite  dimensional,  scalar  transfer  function  for  i  = 
0, 1, . . . ,  n  —  1.  Matrices  of  the  form  given  by  Equation  5.9  are  said  to  be  circulant 
[12],  A  finite  dimensional  system  is  spatially  invariant  if  and  only  if  its  associated 
plant  matrix  is  circulant. 

Define  z  to  be  the  circular  shift  operator,  i.e. 


z  :  Xi  i — y 


xi+1  i^n-1 
x0  i  =  n  —  1 


Note  that  as  a  direct  result  of  this  definition, 


(5.10) 


k  _  k  mod  n 

£  -  A/ 


(5.11) 


for  all  k  G  1L.  Define  the  circular  z  transform  of  the  n-vector  x(s)  to  be 

71—1 

X{s,z)  =  ^2,Xi{s)z\  (5.12) 

i= 0 

Then,  because  of  the  spatial  invariance  property  of  the  plant  matrix  P(s),  the 
system  given  by  Equation  5.8  can  be  written  in  the  z  domain  as 


Y(s,z)  =  G(s,z)U(s,z).  (5.13) 

Here  G(s,z)  is  the  z  transform  of  g(s)  =  [<7o(s)j<7i(s))  •  •  •  )fi,n-i(s)]T- 

The  discrete  Fourier  transform  (DFT)  of  an  n  dimensional  vector  is  equivalent 
to  evaluating  the  z  transform  at  n  evenly  spaced  points  on  the  unit  circle.  In  other 
words,  the  kth  element  of  the  DFT  of  x  is  X (e_j27rfc/n).  As  a  result,  the  DFT  can 
be  used  to  diagonalize  the  system  of  Equation  5.8.  The  resulting  system  is  written 
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as 


i) 

G(s,  1)  0 

0 

U(s,  1) 

Y(s,ei) 

= 

0  G(s,ei) 

0 

U(s,en)) 

V(  ^(n-1  K 

_Y(s,e  n  )  ^ 

0 

0  G(S,e'V’)J 

[  U(s,  eJ(  «  ')  J 

(5.14) 


To  put  this  idea  into  the  context  of  plant  diagonalization  as  discussed  in  the 
previous  chapters,  consider  the  coordinate  transformations  y(s)  =  J-ny(s)  and 
u(s)  =  Pnu(s).  For  a  circulant  plant  matrix  P(s)  we  have 

y(s )  =  P(s)u(s) 

P~ly{s)  =  P(.s)P~lu{s) 

y(s)  =  PnP^p-'uis).  (5.15) 

It  is  not  difficult  to  see  that  the  matrix  PnP(s)P~l  is  equal  to  the  diagonal  matrix 
on  the  right  hand  side  of  Equation  5.14.  Hence,  the  DFT  matrix  can  be  used  to 
diagonalize  any  circulant  plant  matrix. 


5.3  Diagonalizing  Circulant  Plants 

In  the  previous  section  we  showed  that  the  discrete  Fourier  transform  can  be  used 
to  diagonalize  plants  which  exhibit  a  spatial  invariance  property.  Here  we  work 
from  this  result  to  find  ROTs  which  diagonalize  these  plants. 

5.3.1  Background:  The  Discrete  Fourier  Transform 

Let  x  =  [x0,  Xi, . . . ,  xn-i]T  be  a  vector  in  Rn  and  let  X(z)  be  the  circular  z 
transform  of  x  as  defined  in  Equation  5.12.  The  DFT  of  x  is  X(z)  evaluated  at 
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n  evenly  spaced  points  on  the  unit  circle  in  the  complex  plane.  Let  X  denote  the 
DFT  of  x.  The  kth  element  of  the  X  is  X(e_-?27rfc/n),  k  —  0, 1,  2, . . . ,  n  —  1.  Using 
the  definition  of  the  circular  z  transform,  we  get 

n—1 

X(Jfe)  =  X(e~j2nk/n)  =  J2xee~j2nki/n. 

(=0 

Hence  the  DFT  can  be  written  as  a  matrix  whose  kth  row  is 

]_  e-j2irk/n  e~j2ir2k/n  .  .  .  e-j2n(n-l)k/n 

We  will  denote  this  matrix  as  JFn,  i.e.  X  =  Tnx. 

Using  Equation  5.17,  we  see  that  the  (£,  k)t.h  element  of  is  given  by 

n—1 

(5.18) 

m= 0 

If  £  =  k,  then  all  of  the  terms  in  the  sum  on  the  right  hand  side  are  equal  to  1 
and  =  n.  If  i  ^  /c,  then  the  terms  in  the  sum  are  evenly  spaced  points 

on  the  unit  circle.  As  a  result,  they  add  up  to  zero  and  \TTH\(k  =  0.  Hence, 
TnT^  =  nil.  Clearly,  the  matrix  Fn  =  Tnl y/n  is  unitary. 

5.3.2  Gradient  Flows  Diagonalizing  Circulant  Matrices 

For  a  given  set  of  ROT  parameters,  we  consider  the  task  of  Ending  the  ROT  vari¬ 
ables  which  most  closely  diagonalize  any  circulant  plant  matrix.  As  we  have  shown 
in  the  previous  sections,  the  DFT  matrix  accomplishes  the  goal  of  diagonalization 
and  is  unitary  when  properly  scaled.  In  Chapter  4  we  showed  how  to  use  a  gradient 
flow  to  find  the  ROT  variables  which  approximate  the  unitary  SVD  factors  of  an 
arbitrary  complex  matrix.  Here  we  take  a  similar  approach  to  find  an  ROT  which 
diagonalizes  any  circulant  matrix. 


(5.16) 


(5.17) 
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Recall  that  in  the  case  of  diagonalizing  a  symmetric  matrix  Hq,  we  constructed 
a  gradient  flow  to  find  the  0  which  minimized  N  —  0Thfo@  ,  where  N  was  a 
diagonal  matrix  with  the  same  dimension  as  Ha.  Here  instead  of  trying  to  diago¬ 
nalize  a  single  matrix  Hq,  we  are  trying  to  diagonalize  the  whole  set  of  circulant 
matrices.  Towards  this  end,  we  notice  that  the  set  of  n  x  n  circulant  matrices 
forms  an  n  dimensional  vector  space.  As  a  result,  there  exist  a  set  of  basis  vectors 
{Eq,  Ei,  E‘>. . . . ,  En- 1 }  such  that  any  circulant  matrix  A  can  be  written 

n— 1 

A  =  22<*iEi,  (5-19) 

i= 0 

with  at  G  1R  for  i  =  0,l,2,...,n  —  1.  The  task  of  diagonalizing  (or  approximately 
diagonalizing)  any  circulant  matrix  can  then  be  posed  as  the  task  of  simultaneously 
diagonalizing  each  of  the  n  basis  vectors.  Also,  each  E,  is  circulant,  so  the  matrix 
Ni  =  FEiFH  is  diagonal.  Hence,  the  problem  of  diagonalizing  any  circulant  matrix 
can  be  posed  as  the  problem  of  finding  the  0  which  minimizes  the  cost  function 

n- 1  2 

J(0)  =  ki  -  ®HEi®  ~  •  (5-20) 

i= 0 

As  in  Chapter  4,  simple  matrix  manipulation  shows  that  minimizing  J  is  equiv¬ 
alent  to  maximizing  the  “diagonalness”  function 

n—1 

$  (©)  =  J2  retr  (A rtH<dHEi<d^  .  (5.21) 

i= 0 

The  unitary  operator  0  is  an  ROT,  i.e. 

L 

0  =  n*0*-  (5.22) 

k= 1 
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Each  ©fc,  k  —  1, 2, . . . ,  L  is  of  the  form 


e\  0  •  •  •  0 

o  el  ■■■  o 

©fc  =  ,  (5.23) 

.0  ■  ■  ■  0  elk  _ 

where  01-  £  U (nkj), 

mk 

^^nkj  =  n  Wk  =  1, 2, . . . ,  L,  (5-24) 

3  = 1 

and  the  Os  represent  appropriately  dimensioned  blocks  of  zeros.  We  define  Mk  = 
U(nk i)  x  U(nk2)  x  •  •  •  x  U{nkmk).  Then  each  ©fc  belongs  the  smooth  Lie  subgroup 
Mk  £  U(n).  The  Lie  algebra  of  Mk  is 

TeMk  =  u(nk i)  ©  u(nfc2)  ©  •  •  •  ©  u(nfemJ.  (5.25) 


A  vector  in  TeMk  can  then  be  written 

c o\  0  •  •  •  0 

0  uf  ■■■  0 

Vlk  =  .  (5.26) 

_  0  ...  0  _ 

where  £  u(nkj )  for  j  —  1,2, ,  mk. 

Here,  we  assume  that  the  configuration  parameters  of  the  ROT  are  given  and 
fixed.  In  other  words,  the  level  L  is  fixed  and  for  k  =  1,2 ,...L,  the  following 
quantities  are  fixed:  Pk,  mk,  and  nkj,  j  =  1,2, .. .  ,mk.  Our  task  is  to  End  the 
ROT  variables,  i.e.  the  O^s,  which  maximize  (j). 


Let  M  =  Mi  x  M2  x  •  •  •  x  ML.  The  set  M  is  a  Lie  group.  An  element  of  M 
is  written  as  the  L-tuple  T  =  (01;  02, . . . ,  O L).  A  vector  in  T^M  is  written  as 
X  =  (©ifli,  02fl2, . . . ,  O l^l),  where  flk  £  TeMk  for  each  k. 
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The  diagonalness  function  (j)  can  now  be  thought  of  as  a  functional  on  M, 
(j)  :  M  — »  H.  As  in  Chapter  4,  our  approach  to  searching  for  the  T  which  maximizes 
(j)  is  to  flow  along  the  gradient  vector  field  of  0  on  M. 

Recall  that  V0(T)  on  M  is  defined  by  the  following  properties: 

1.  V0(T)  e  T^M  VT  e  M. 

2.  D^(X)  =  (V0(T),  X)  VX  G  T*M. 

The  hrst  property  states  that  VbfT)  in  contained  in  the  tangent  space  of  M  at 
the  point  \F  This  means  that  V0(T)  must  be  of  the  form 

V0(T)  =  (e^,  02o^, . . . ,  0Ln^  ,  (5.27) 

where  G  TeMk,  k  =  1,2,  For  convenience  of  notation,  we  have  sup¬ 

pressed  the  fact  that  Oj0  depends  on  \F 
From  the  second  property,  we  must  have 

. . . ,  ©l^l))  =  (V0(©i, . . . ,  ©L),  (©ifii, . . . ,  ©l^l)) 

=  ((0!^, . . . ,  ©Ln^),  (©!«!, . . . ,  eLnL)) 

=  J]retr 

fc=i 

L 

=  —  retr  (5.28) 

fc=i 

Finding  an  expression  for  D(j)y(X): 

D<h(X)  = 

retr  (eC  EEo  ((ntl  fie<)  (rt*+1  e,et))" s.  at.  ?A) 

+^f  cn?=i  (nfc,1  p&t)  ns*  (n?_s+1  ^e*))  • 

(5.29) 
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To  simplify  notation,  we  define  two  sequences  of  recursive  orthogonal  trans¬ 


forms: 


(5.30) 


e„  =  n  p,e> 


(5.31) 


Using  this  notation,  we  get 


D<fa(X)  =  retrl  E&  +  N? ®TEiQknk® 


k= 1  i= 0 


retr  E  E  NfOMVkES  +  NlHQHEiQknket 


k= 1  i= 0 


E  retr  E  Ei®k  1  fi 


(5.32) 


Claim  5.3.1  Let  {_E0,  £/,  . . . ,  £^n_i}  6e  an  orthonormal  basis  for  the  space  ofnxn 
circulant  matrices.  Let  Ni  =  FEiFH ,  /or  z  =  0, 1, . . . ,  n  —  1,  where  F  is  the  DFT 
matrix  scaled  by  1/ sjn.  Then  the  matrix 


r*  =  £[6‘JV"6»'e*£;‘e, 


(5.33) 


is  shew-hermitian  for  k  =  1,2 , ...  ,L. 


Proof: 


Expanding  the  Lie  bracket  and  substituting  for  Ni:  we  get 


Tfc  =  Y,®KFHE?F®kQkEiQk-QkEiQk®kFHE?F®\ 


QkFH  Y,E?F'9%e%EiQk9kFI 


-F®«®»El®k®kFHE?)F®«. 


(5.34) 
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Define  A  =  Ou®kFH  and  define 


n— 1 


K  =  Yl  EfFQ^e^EiQkQkF11  -  F0% ©f  EtQkQkFH E^1 
2=0 
n— 1 

=  E?AHEiA  -  AuEiAEf.  (5.35) 

i= 0 

Substituting  K  into  Equation  5.34  and  taking  the  (Hermitian)  transpose  gives 


r f  =  (  QuFM  KFQ 


H  TS  paH 


=  ®kFHKHF®?. 


(5.36) 


As  a  result,  1^.  is  skew  if  and  only  if  K  is  skew.  By  definition,  K  is  skew  if 
K  +  Kh  =  0.  Finding  an  expression  for  K  +  KH: 

n— 1  / n— 1  y  H 

K  +  Kh  =  Y  FfAHEiA  -  AHEiAEf  +  I  Y  E?  AH  EtA  -  AHEiAEf 


i= 0 
n—1 


.  i=0 


Y  E?AHEiA  -  AHEiAE f  +  AH  E?  AE.t  -  E.tAH E? A.  (5.37) 


t=0 


Since  { E,}  is  a  basis  for  circulant  matrices  and  Ef  is  a  circulant  matrix,  we 
can  write  each  Ef1  as 

n—1 

Ei  =Y^aijEj,  (5.38) 

3=0 

where  each  u%3  G  1R.  Since  { Et }  is  orthonormal,  we  have  an  explicit  expression  for 


Otij\ 


aij  ~  {Ei  >Ej} 

=  retr  (EiEj) . 


(5.39) 


Note  that  retr  {EiEj)  =  retr  (. EjEi ),  so  atJ  =  ay,  for  each  i  and  j.  Substituting  for 
Ef  in  Equation  5.37  yields 
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K  +  Kh 

n— 1  / n— 1  \  / n— 1 

=  E  E^£dA"B‘A-A"B‘A  E“«B.f 

*=0  \j=0  )  \j= o 

AEi  —  EiAH  ^  °i] Ej ^  A 

n— 1  n— 1 

=  EE««  ((E]AhE1A  -  EiAHEjA )  +  (AHEjAEi  -  AHEiAEj)) 

i= 0  j=0 
n— 1  n— 1 

=  E  E  “«  {{(E3AHEtA  -  EiAHEjA)  +  (. AHEjAEi  -  AHEiAEj)) 

i— 0  j=i-\- 1 

+cr/i  ((EiAHEjA  -  EjAHEiA )  +  [AHEiAEj  -  AHEJAEl))) 

n— 1 

+  a*  [[EiAH Ei A  -  E1AhE,A)  +  (AH  EiAEt  -  AH EtAE,)) 

i= 0 

n— 1  n— 1 

=  E  E  (“« -  aji)  ((EjAHEiA  -  EiAHEjA) 

i— 0  j=i+l 

+  ( AHEjAEi  -  AHEiAEj))  .  (5.40) 

But  we  have  already  shown  that  a,,  =  am ,  so  K  +  A  ^  =  0.  ■ 


Looking  back  at  Equation  5.28,  we  see  that  in  order  to  have  a  valid  gradient 
flow,  each  Qk  in  Equation  5.32  must  be  multiplied  on  the  left  by  — where 
^v</>  g  Te(Mj.).  Claim  5.3.1  tells  us  that  the  matrices  which  multiply  the  Qj.s 
on  the  left  are  skew-hermitian.  To  put  them  in  Te[Mk)  we  introduce  a  set  of 
natural  projection  operators,  II/,.  :  TeU(n )  — *  Te(Mk).  Lb,  projects  from  the  set  of 
skew-symmetric  matrices  to  the  set  of  block  diagonal  skew-symmetric  matrices  by 
setting  all  off  diagonal  blocks  to  zero.  Now  we  can  use  Fact  4.1.1  to  state 

L  /  /  n— 1 


DMX)  =  5>etr  nfc  ^  ©fcATfef,©^© 


A  i 


k= 1 


,  i=0 


Now 


'  n—  1 


J= 0 


Bit  e.ivf  ef ,  ©f EiQk 


e  TeMk 


(5.41) 


(5,42) 
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for  each  k,  so  the  gradient  flow  which  satisfies  both  properties  is  given  by  the 


assignment 


'  n— 1 


*%+  =  -nfc  £  n* 


,  i=0 


QkN*B*,Q*EiQ 


The  gradient  flow  is  then  given  by  the  couple  matrix  ODEs 


'71—1 


Qk  =  -oknk  (  ^  1 0fciyf0f,0fE80A 

for  k  —  1, 2, . . . ,  L. 


,  ©fc(0)  —  0fcO, 


(5,43) 


(5.44) 


„  i=0 


5.3.3  Example 

Here  we  apply  this  plant  diagonalization  technique  to  the  example  of  a  flexible 
beam  with  periodic  boundary  conditions.  This  system  can  be  thought  of  as  a 
flexible  ring.  The  model  we  use  is  the  beam  model  described  and  modeled  in 
Appendix  B  with  the  boundary  conditions  at  the  ends  of  the  beam  replaced  by 
periodic  boundary  conditions.  In  this  example,  we  consider  a  beam  with  8  inputs 
and  8  outputs.  The  sensors  and  actuators  are  evenly  spaced  around  the  ring  so 
that  the  system  has  the  spatial  invariance  property. 

The  plots  in  Figure  5.1  shows  the  responses  of  the  8  outputs  to  an  impulse 
applied  to  the  first  input.  Figure  5.2  depicts  the  circulant  nature  of  the  plant.  In 
this  figure,  the  intensity  of  the  (i,  j)th  pixel  represents  the  energy  contained  in  the 
response  of  the  ith  output  to  an  impulse  applied  to  the  j th  input.  We  call  this 
matrix  the  impulse  energy  matrix. 

We  found  diagonalizing  basis  transformations  as  described  in  this  chapter  using 
3,  5,  and  7  level  ROTs.  Figures  5.3,  5.4,  and  5.5  show  the  impulse  energy  matrices 
for  these  transforms.  The  coordinate  transforms  are  complex  so  the  resulting 
impulse  responses  are  complex  valued.  To  compute  the  impulse  energy,  we  take 
the  magnitude  of  the  impulse  response. 
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Figure  5.2:  Impulse  energy  matrix  for  untransformed  plant. 
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Impulse  Energy  Matrix  of  Flexible  Ring  with  Level  3  ROT 


Figure  5.3:  Impulse  energy  matrix  for  plant  transformed  with  level  3  ROT. 


Impulse  Energy  Matrix  of  Flexible  Ring  with  Level  5  ROT 


Figure  5.4:  Impulse  energy  matrix  for  plant  transformed  with  level  5  ROT. 
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Chapter  6 


Conclusions 


In  this  thesis,  we  have  presented  the  foundations  of  a  novel  approach  to  signal 
processing  and  control  on  large  distributed  control  networks.  Onr  approach  is 
to  apply  efficient  data  transformations  to  the  input  and  output  vector  of  a  plant 
in  an  effort  to  make  the  resulting  plant  matrix  look  diagonal.  In  Chapters  2 
and  4  we  proposed  efficient  transforms  capable  of  “approximately”  diagonalizing 
any  constant  matrix.  These  ideas  are  useful  for  stable  linear  systems  in  which 
the  transient  behavior  can  be  ignored  as  well  as  systems  which  exhibit  strong, 
sharp  peaks  (resonances)  in  their  frequency  response.  In  Chapter  5  we  developed 
a  method  of  diagonalizing  dynamic  plants  which  exhibit  a  linear  spatial  invariance 
property. 

We  conclude  the  thesis  with  a  quantitative  comparison  of  the  matrix  diagonal- 
ization  methods  presented  followed  by  some  suggestions  for  future  research. 


121 


6.1  Comparison  of  Proposed  Approximate  Diag- 
onalization  Methods 

In  this  section  we  compare  and  contrast  the  methods  of  approximate  matrix  di- 
agonalization  that  we  have  developed  in  this  thesis.  Onr  purpose  is  not  only  to 
highlight  the  merits  of  the  proposed  methods,  bnt  also  to  point  out  where  there  is 
room  for  improvement. 

6.1.1  Performance  Measures 

Here  we  define  the  performance  measures  that  will  be  used  to  compare  the  var¬ 
ious  data  transformation  techniques.  The  costs  we  consider  are  communication 
and  computational  costs  associated  with  transforming  the  output  vector.  We  also 
define  a  performance  measure  which  indicates  how  well  the  proposed  transforms 
diagonalize  the  plant  matrix. 

Global  Data  Transfers 

One  common  network  configuration  is  the  common  bus  topology.  In  this  configu¬ 
ration,  all  of  the  nodes  on  the  network  share  one  serial  communication  bus.  The 
cost  of  a  data  transfer  in  this  configuration  does  not  depend  on  the  positions  of 
the  sending  receiving  nodes.  It  is  as  expensive  for  a  node  to  send  information  to 
its  neighbor  as  it  is  to  send  information  to  a  node  on  the  other  side  of  the  network. 
In  fact,  for  the  same  cost,  the  sending  node  can  communicate  its  data  to  every 
node  on  the  network.  We  call  any  such  communication  a  global  data  transfer. 

As  an  example,  consider  a  1-level  ROT  on  a  linear  array  as  described  in  Sec¬ 
tion  3.2.1.  Recall  that  the  first  (and  only)  communication  step  is  for  each  evenly 
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indexed  node  to  send  the  value  of  its  sensor  output  to  the  node  to  its  left.  For  an 
array  with  n  nodes,  this  requires  n/2  global  data  transfers. 

Nearest  Neighbor  Data  Transfers 

To  measure  this  communication  cost,  we  assume  that  each  node  on  the  network  can 
exchange  information  only  with  its  nearest  neighbors.  The  transfer  of  data  from 
a  node  to  its  neighbor  is  called  a  nearest  neighbor  data  transfer.  Multiple  data 
transfers  are  required  to  exchange  information  between  non-neighboring  nodes. 
For  example,  in  order  to  send  a  piece  of  data  from  the  ith  node  of  a  linear  array 
to  the  j th  node,  a  total  of  |i  —  j\  nearest  neighbor  data  transfers  are  required. 

Nearest  Neighbor  Transfer  Time 

In  a  nearest  neighbor  network,  independent  data  transfers  can  be  performed  at  the 
same  time.  We  define  the  nearest  neighbor  communication  time  to  be  the  total 
time  required  to  carry  out  all  communications,  taking  into  account  the  fact  that 
many  can  be  carried  out  simultaneously.  The  unit  used  to  measure  this  time  is  the 
amount  of  time  required  to  carry  out  one  nearest  neighbor  data  transfer.  We  call 
this  time  a  step. 

To  send  a  piece  of  data  from  the  ith  node  of  a  linear  array  to  the  j  th  node,  a 
total  of  i  —  j  nearest  neighbor  data  transfers  are  required.  These  transfers  cannot 
be  performed  in  parallel,  so  the  communication  time  is  |i  —  j\  steps.  On  the  other 
hand,  a  1-level  ROT  performed  on  a  linear  array  with  n  nodes  requires  n/2  nearest 
neighbor  communications.  All  of  these  can  be  performed  simultaneously,  so  the 
corresponding  communication  time  is  1  step. 
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Implementation  Multiplications 


Here,  we  measure  computational  complexity  by  counting  the  total  number  of  mul¬ 
tiplications  required  to  implement  the  output  vector  transformation. 

Implementation  Multiplication  Time 

Here  we  measure  the  amount  of  time  required  to  carry  out  the  multiplications  for 
the  output  vector  transformation.  In  this  measure,  we  take  into  account  the  fact 
that  some  operations  can  be  performed  simultaneously.  The  unit  of  measurement 
is  the  time  required  to  complete  one  multiplication.  We  call  this  time  a  step. 

Off-line  Multiplications 

In  addition  to  implementation  costs,  there  is  a  cost  associated  with  finding  the 
transformation  to  be  applied.  We  measure  this  cost  by  counting  the  number  of 
multiplications  required.  We  do  not  put  as  much  weight  on  this  cost  as  we  do  on 
the  implementation  costs;  the  philosophy  of  this  thesis  is  that  it  is  worth  the  extra 
off-line  effort  to  improve  the  on-line  costs.  We  include  this  cost  in  our  comparison 
to  give  an  indication  of  just  how  much  extra  off-line  work  is  required. 

Diagonal  Energy  Percentage 

To  help  compare  the  performances  of  the  various  methods  presented  in  this  thesis, 
we  develop  a  measure  that  indicates  the  “diagonalness”  of  an  approximately  diag¬ 
onalized  matrix.  Let  H  be  an  n  x  m  matrix.  Let  x  be  the  vector  whose  elements 
are  the  diagonal  elements  of  H ,  i.e.  x(i)  =  [H]a  for  %  =  1,2, ... ,  min(n,  m).  We 
define  the  diagonal  energy  percentage  to  be 

v( H)  =  S  x  100,  (6.1) 
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where  ||if||  denotes  the  Frobenius  norm  of  H. 

6.1.2  Implementation  of  Exact  SVD  Factors 

In  Section  2.1  we  pointed  out  that  the  SVD  factors  can  be  used  to  diagonalize 
a  real-valued  matrix.  We  also  claimed  that  they  require  too  much  computation 
and  communication  to  be  implemented  in  real  time  on  a  large  distributed  control 
network.  Here  we  briefly  discuss  such  an  implementation  of  the  SVD  factors  so 
that  we  can  quantify  these  remarks  and  compare  the  performance  of  SVD  factors 
to  the  other  matrix  diagonalization  techniques  proposed  in  this  thesis.  We  discuss 
only  the  implementation  of  the  output  vector  transformation.  The  implementation 
of  the  input  vector  transformation  can  be  inferred  from  this  discussion. 

The  most  straight  forward  method  is  what  we  refer  to  as  the  centralized  im¬ 
plementation.  In  this  implementation,  each  network  node  sends  its  piece  of  the 
output  vector  to  a  central  location  where  the  computations  are  carried  out.  If  we 
assume  that  the  n  nodes  on  the  network  share  a  common  bus  allowing  them  to 
communicate  with  every  node  on  the  the  network,  then  the  communication  step 
requires  (n  —  1)  data  transfers.  If  we  instead  assume  that  only  nearest  neighbor 
communication  is  required,  then  the  fewest  data  transfers  required  to  get  all  of  the 
data  to  a  central  location  is 

n/ 2-1 

Nc  =  2  k.  (6.2) 

k= 1 

Once  all  of  the  data  is  collected,  the  computation  which  needs  to  be  executed 
consists  of  the  multiplication  of  an  n  x  n  matrix  with  an  n  dimensional  vector. 
This  requires  n 2  multiplications  and  n(n  —  1)  additions. 

The  SVD  factors  can  also  be  realized  using  a  parallel  implementation.  Here 
each  node  on  the  network  must  send  its  piece  of  the  data  vector  to  every  other  node 
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on  the  network.  Then  the  matrix  multiplication  is  carried  out  in  parallel,  with  each 
node  computing  the  multiplication  of  one  row  of  the  matrix  with  the  data  vector. 
If  we  assume  that  the  n  nodes  on  the  network  share  a  common  bus  allowing  them 
to  communicate  with  every  node  on  the  the  network,  then  the  communication  step 
requires  (n  —  1)  data  transfers.  If  only  nearest  neighbor  communication  is  allowed, 
then 

71—1 

Np  =  2j2k  (6.3) 

k= 1 

data  transfers  are  required.  Once  all  of  the  necessary  transfers  have  been  made, 
each  node  performs  the  multiplication  of  one  row  of  the  SVD  factor  by  the  data 
vector.  The  total  number  of  operations  necessary  to  complete  the  transform  is  the 
same  as  in  the  centralized  implementation,  but  since  the  parallel  implementation 
allows  the  row  multiplications  to  be  carried  out  simultaneously  the  computation 
time  is  reduced.  Hence,  the  computations  are  completed  in  the  time  it  takes  to 
perform  n  multiplications  and  (n  —  1)  additions. 

6.1.3  Comparison  Tables 

Table  6.1.3  lists  these  measures  for  a  variety  of  approximate  matrix  diagonalization 
techniques  applied  to  the  plant  matrix  of  a  flexible  membrane  driven  by  a  linear 
array.  Table  6.1.3  compares  the  same  techniques  applied  to  the  plant  matrix  of  a 
flexible  beam  at  its  sixth  resonance.  In  these  tables,  the  SVD  measures  correspond 
to  the  parallel  implementation  of  the  SVD  factors. 

Though  the  costs  and  performances  of  the  various  transforms  vary  from  exam¬ 
ple  to  example,  we  can  make  some  generalizations  about  the  comparisons.  The 
wavelet  packet  transform  generally  require  the  lowest  number  of  implementation 
multiplications,  the  lowest  multiplication  time,  and  the  lowest  off-line  multiplica- 
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SVD 


Wavelet 


ROT  10  ROT  20 


ROT  31 


ii 


Global  Data 

Transfers 

32 

110 

160 

320 

496 

Nearest  Neighbor 

Data  Transfers 

992 

480 

160 

320 

496 

Nearest  Neighbor 

Transfer  Time 

62  steps 

40  steps 

10  steps 

20  steps 

30  steps 

Implementation 

Multiplications 

1024 

252 

640 

1280 

1984 

Implementation 

Mult.  Time 

32  steps 

20  steps 

40  steps 

80  steps 

124  steps 

Off-line 

Multiplications 

~  3  x  104 

~  6  x  103 

-  1010 

~  1011 

~  1014 

Diagonal 

Energy  % 

100% 

93.4% 

64.2% 

88.7% 

99.9% 

Table  6.1:  Comparison  table  for  membrane  example. 
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SVD 


Wavelet 


ROT  10  ROT  20 


ROT  31 


ii 


Global  Data 

Transfers 

32 

104 

160 

320 

496 

Nearest  Neighbor 

Data  Transfers 

992 

218 

160 

320 

496 

Nearest  Neighbor 

Transfer  Time 

62  steps 

26  steps 

10  steps 

20  steps 

30  steps 

Implementation 

Multiplications 

1024 

240 

640 

1280 

1984 

Implementation 

Mult.  Time 

32  steps 

20  steps 

40  steps 

80  steps 

124  steps 

Off-line 

Multiplications 

~  3  x  104 

~  6  x  103 

-  1010 

~  1011 

~  1014 

Diagonal 

Energy  % 

100% 

60.6% 

54.3% 

69.7% 

99.3% 

Table  6.2:  Comparison  table  for  the  example  of  a  flexible  beam  at  its  6th  resonance. 
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tions.  The  ROTs  exhibit  the  lowest  nearest  neighbor  transfer  time.  The  higher 
the  level  of  the  ROT,  the  better  its  performance  in  terms  of  diagonal  energy,  but 
this  performance  comes  at  the  cost  of  every  other  measure  listed. 

6.2  Suggestions  for  Future  Work 

The  main  theoretical  question  that  remains  is  the  issue  of  the  convergence  of  the 
ROT  gradient  flow  equations  and  the  characteristics  of  the  equilibrium  points 
to  which  these  equations  converge.  We  begin  to  address  these  issues  in  for  the 
case  of  transforming  a  4  dimensional  data  vector  (See  Appendix  A),  but  a  more 
thorough  understanding  is  required.  Such  an  understanding  may  yield  bounds  on 
the  performance  of  an  ROT  or  algorithms  to  select  the  “best”  set  of  permutation 
matrices. 

Another  issue  that  needs  to  be  addressed  is  the  numerical  integration  of  the 
ROT  gradient  flow  equations.  The  computations  required  to  find  a  31-level,  32  x  32 
ROT  are  almost  prohibitive.  More  efficient  numerical  techniques  are  necessary  if  we 
are  to  find  ROTs  for  larger  systems.  This  boils  down  to  the  problem  of  integrating 
the  nonlinear  matrix  ODE 

0  =  00(0),  (6.4) 

on  0(n)  where  0(0)  G  o(n)  for  every  0  G  0(n).  It  is  important  to  use  a  numerical 
integration  scheme  which  respects  the  constraint  that  0  remains  orthogonal.  The 
recursive  sequence 

0fc+i  =  0fcexp(aO(0))  (6.5) 

provides  one  such  integration  method.  This  method  is  analogous  to  the  Euler 
method,  where  the  scalar  a  is  the  step  size.  Brockett  [7]  and  Moore,  Mahony,  and 
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Helmke  [29]  have  developed  recursive  methods  for  their  work  which  combine  this 
idea  with  a  clever  choice  of  a.  These  methods  give  fast  convergence  to  the  desired 
solution.  A  similar  approach  should  be  possible  for  the  ROT  equations.  Other  nu¬ 
merical  gains  might  be  found  by  using  generalizations  of  the  Cayley  transform  (see 
Tsiotras,  Junkins,  and  Schaub  [39])  to  approximate  the  computationally  intensive 
matrix  exponential. 

Throughout  the  thesis,  we  have  assumed  that  the  ROT  parameters  are  fixed 
and  given.  In  particular,  we  selected  the  parameters  based  on  a  nearest  neighbor 
communication  scheme.  Almost  certainly,  better  performance  can  be  achieved  by 
changing  the  ROT  parameters.  One  idea  towards  this  end  is  to  use  devise  permu¬ 
tation  matrices  which  mimic  the  communication  schemes  of  the  transforms  in  the 
wavelet  packet.  Alternatively,  it  may  be  possible  to  use  more  sophisticated  search¬ 
ing  algorithms  to  optimize  the  ROT  variables  and  parameters  simultaneously. 
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Appendix  A 


Analysis  of  ROTs  on  SO  (A) 


Here  we  study  the  some  of  the  properties  of  a  recursive  orthogonal  transform 
(ROT)  composed  of  4  x  4  matrices.  The  aim  of  this  exercise  is  to  provide  some 
insight  into  the  nature  of  the  ROT  and  perhaps  lay  some  groundwork  for  future 
research.  Specifically,  we  examine  the  3-Level  ROT 

O  =  OiP2©2-P303-  (A.l) 


Here, 


0\  02x2 

02x2  $2 

where  6 *  G  SO (2)  for  i  <G  {1,2,3},  j  G  {1,2}.  The  permutation 


(A.2) 


matrices  are 


0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

1 

0 

1 

0 

0 

0 

p2  = 

0 

0 

0 

1 

and  P3  = 

0 

1 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

(A.3) 


Equation  A.l  can  be  implemented  as  a  data  transformation  on  a  four  node  linear 
array  using  only  nearest  neighbor  communication  (See  Section  3.2.1).  Matrices  of 
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the  form  given  in  Equation  A. 2  form  a  compact  subgroup  (torus)  of  AO(4).  We 
will  denote  this  subgroup  by  M. 

Each  component  in  the  product  on  the  right  hand  side  of  Equation  A.l  is  in 
SO  (A),  so  0  is  also  in  AO(  4).  The  question  that  naturally  arises  at  this  point  is 
the  following:  what  subset  of  AO (4)  can  be  represented  by  the  above  ROT?  Or, 
more  to  the  point,  is  it  possible  to  represent  any  $  G  AO (4)  with  a  product  of  the 
form  given  in  Equation  A.l? 

We  know  that  the  dimension  of  AO  (4)  is  six.  In  order  to  have  any  chance  of 
representing  all  of  AO  (4),  0  must  have  at  least  six  degrees  of  freedom  which  are, 
in  some  sense,  independent.  Each  0,  in  Equation  A.l  has  two  degrees  of  freedom, 
so  O  has  a  total  of  six.  It  is  not  clear  that  these  six  dimensions  are  independent, 
though.  For  example,  if  P2  and  P:>  were  both  the  identity  matrix,  then  0  could 
only  represent  elements  in  the  two  dimensional  subgroup  M  C  AO(4). 

To  answer  these  questions  we  look  to  the  theory  of  symmetric  subalgebras. 
Here  we  give  a  brief  presentation  of  some  of  these  tools  and  then  show  how  to 
apply  them  to  the  above  ROT. 

A.l  The  KAK  Representation 

Here  we  present  a  brief  introduction  to  the  KAK  representation  of  semisimple  Lie 
groups.  We  follow  Hermann  [23]  and  present  the  main  theorems  without  proof. 
Readers  interested  in  the  details  of  this  material  are  referred  to  Hermann  or  Hel- 
gason  [20]. 

Definition  A. 1.1  Let  G  be  a  connected,  semisimple  Lie  group  with  finite  center. 
Let  g  be  the  Lie  algebra  ofG.  Let  £  be  a  subspace  of  g  and  let  p  be  the  complement  of 
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1  in  g.  Then  f  is  said  to  be  a  symmetric  subalgebra  of  q  if  the  following  bracketing 
conditions  hold: 


1.  [M]  C  t 

2.  [t,p]cp. 

3.  [p,p]  C  t 


Associated  with  the  subalgebra  1  C  g  is  a  subgroup  K  C  G.  The  subgroup  K 
is  said  to  be  a  symmetric  subgroup  of  G.  Since  G  has  a  finite  center,  K  is  compact. 

Since  K  is  compact,  the  exponential  map  maps  f  onto  K  [33].  As  a  result,  for 
any  k  G  K  there  exists  X  £  t  such  that  k  =exp(X). 

Theorem  A. 1.1  (Theorem  6—4  from  Hermann  [23]) 

Let  G,  K ,  g,  f,  and  p  be  as  defined  above.  Let  a  be  a  maximal  Abelian  subalgebra 
of  p  and  let  a'  be  any  Abelian  subalgebra  ofp.  Then  there  exists  k  G  K  such  that 

ka'k -1  C  a.  (A. 4) 

Theorem  A. 1.2  (Theorem  6—5  from  Hermann  [23]) 

Let  G,  K,  0,  t,  and  p  be  as  defined  above.  Let  P  be  the  image  of  p  under  the 
exponential  map.  Then  any  g  G  G  can  be  represented  as 

g  =  pk,  (A. 5) 


where  p  G  P  and  k  G  K. 

These  two  theorems  are  combined  to  give  us  the  so-called  KAK  representation 
of  the  group  G.  First,  Theorem  A. 1.2  says  that  any  g  G  G  can  be  represented  as 
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exp (Xp)k  where  Xp  G  p  and  k  G  K.  Now  let  a  be  a  maximal  Abelian  subalgebra 
of  p.  A  direct  result  of  Theorem  A. 1.1  is  that  every  element  of  p  can  be  written 
Xp  =  kiXak^ 1  where  k\  G  K  and  Xa  G  a.  Substituting,  we  get 

g  =  exp(k1Xak41)k 

=  kiexp(Xa)k41k.  (A. 6) 

Since  K  is  a  group,  k f 1  k  G  K.  Letting  A  denote  the  image  of  a  under  the 
exponential,  we  see  that  any  g  G  G  can  be  written  as 

g  =  k4ak2,  (A. 7) 

where  k { .  k2  G  K  and  a  G  A.  This  is  known  as  the  KAI\  representation. 

A. 2  ROT  Representation  of  50(4) 

To  begin,  we  note  that  <50(4)  is  a  simple  group  [37],  which  means  it  is  also  semisim¬ 
ple.  As  a  result,  the  theory  developed  in  the  previous  section  can  be  applied.  Next, 
we  define  the  following  basis  vectors  for  50  (4): 


0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

0 

Ai  — 

1 

0 

0 

0 

A2  — 

0 

0 

-1 

0 

A3  — 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

A4  = 

0 

0 

0 

0 

A-5  = 

0 

0 

0 

-1 

A-6  = 

0 

0 

0 

0 

0 

0 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

1 

0 

0 

1 

0 

0 

0 

(A.8) 
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Since  M  is  compact,  any  ©j  G  M  can  be  written  as  the  exponential  of  an  element 
in  the  Lie  algebra  of  M.  Using  the  above  basis  vectors,  we  can  write 

©i  =  exp(oqAi  +  ol  2A4),  (A. 9) 

where  a\  and  a2  are  real,  i  G  {1,2,3}.  Substituting  into  Equation  A.l,  0  can  be 
written 

0  =  exp(oqAi  +  a\Ai)P2ey^(ya2lAi  +  alA4)P3exp(al  Ai  +  alA4) 

=  exp(a}A.!  +  a\A4)P2ex\i(KalAi  +  a\A4)P2  eyi^a^Ai  +  a\A4) 

=  exp(ajA.!  +  alA4)exp(alP2AiP2  +  a\P2A4P2)ex\){a\Ai  +  a|Ai) 

=  exp(a}A.!  +  «2A4)exp(— a\A§  +  «2A2)exp(a'j!A1  +  a\A4).  (A. 10) 

Let  1  =  span{Ai,A4}.  It  is  easy  to  check  the  bracket  conditions  in  Defini¬ 
tion  A.  1.1  to  see  that  Ms  a  symmetric  subalgebra  of  so  (4).  Let  p  be  the  complement 
of  t  in  so(4),  i.e.  p  =  span {A2,  A3,  A5,  Aq}.  Explicitly,  the  bracket  [Ai,A4]  =  0 
implies  that  [6,  6]  C  t.  The  brackets 


[Ai,  A2]  =  —A3 

[A4,  a2] 

—  Aq 

[Ai,  A3]  =  A2 

[A4,  a3] 

=  a6 

(A. 11) 

[Ai-,  A5]  =  —Aq 

[A4,  Aq] 

=  — A2 

[Al,  Aq]  =  Aq 

[A4,  Aq] 

=  — a3 

imply  that  [fi,p]  C  p.  Finally,  the  brackets 

[A-2,  A3]  =  —  Ai 

[A3,  Aq] 

=  0 

[A2,  A5]  =  —  A4 

[A3,  Aq] 

=  A4 

(A- 12) 

[A2,  Aq]  =  0 

[Aq,  Aq] 

=  — Ai 

imply  that  [p,  p]  C  t  From  the  last  set  of  brackets,  is 

also  clear  that  the  subalgebra 

a  =  span{A2,  Aq}  is  a  maximal  Abelian  subalgebra  of  p. 
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Using  the  I\AK  representation,  we  see  that  any  $  G  SO( 4)  can  be  written  as 


$  =  exp(Xfci)exp(Xa)exp(Xfc2),  (A.13) 

where  Xki,Xf-2  £  f  and  Xa  G  a.  This  expression  is  of  exactly  the  same  form  as 
the  expression  for  0  given  in  Equation  A.  10.  Therefore,  the  ROT  O  can  be  used 
to  represent  any  $  G  SO  (A). 

This  representation  exploits  a  somewhat  unique  feature  of  so  (4):  its  contains 
Abelian  subalgebras  which  are  also  symmetric.  This  also  holds  for  50 (3),  but  it 
does  not  in  general  hold  for  so(n)  where  n  >  4.  This  means  that  the  above  analysis 
does  not  easily  extend  to  address  ROTs  on  higher  dimensional  spaces. 

A. 3  Convergence  of  ROT  Equations 

In  the  previous  sections  we  showed  that  a  3-level  ROT  given  in  Equation  A.l  can 
be  used  to  represent  any  matrix  in  S'0(4).  This  means  that  given  any  symmetric 
4x4  matrix  Hq,  there  exists  a  0  of  the  form  given  in  Equation  A.l  so  that  the 
matrix  H  =  077/()0  is  diagonal.  However,  this  does  not  necessarily  mean  that 
the  gradient  flow  equations  we  derived  to  End  the  ROT  variables  in  Chapter  4 
converge  to  O.  Hence,  it  is  necessary  to  study  the  convergence  properties  of  the 
ROT  gradient  flow  equations. 

As  we  pointed  out  in  Chapter  4,  the  system  of  ROT  equations  evolves  on  a 
compact  manifold.  By  the  properties  of  a  gradient  flow  on  a  compact  manifold,  the 
solution  to  the  ROT  equations  exists  for  all  time  and  converges  to  an  equilibrium 
point  of  the  system.  The  first  step  is  understanding  the  convergence  properties  of 
the  system  of  ROT  equations  is  to  investigate  the  nature  of  the  equilibrium  points 
of  the  system.  Unfortunately  this  is  a  tedious  task,  and  even  in  the  relatively 
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simple  case  of  4  x  4  matrices  it  is  difficult  to  get  clear  answers.  Here  we  analyze 
the  equilibrium  points  of  the  ROT  equations  on  4  x  4  matrices  for  1 — ,2 — ,  and 
3-level  ROTs.  In  particular,  we  study  the  characteristics  of  the  approximately 
diagonalized  matrix  which  result. 

A. 3.1  1-Level  ROT 

For  a  1-level  ROT,  we  simply  have 

0  =  ©!,  (A. 14) 

where  ©i  is  of  the  form  given  in  Equation  A. 2.  Following  Section  4.1.2,  the  gradient 
flow  equation  to  search  for  the  ROT  variables  which  most  nearly  diagonalize  the 
symmetric  matrix  Ho  is 

0!  =  -©in  [TV,  ©ftf0©i]  ,  0i(O)  =  010,  (A. 15) 

where  N  is  a  4  x  4  diagonal  matrix  with  distinct  entries  and  n  is  the  projection 
operator, 

On  012  Ul3  a14  all  a12  0  0 

<221  °22  <223  a24  <221  a22  0  0 

<231  ^32  <233  ®34  0  0  (Z33  <234 

&41  $42  &43  &44  J  0  0  (Z43  (Z44 

Let  ©1*  be  an  equilibrium  point  of  Equation  A. 15  and  let  H  =  ©^i70©i*  be  the 
resulting  approximately  diagonalized  matrix.  Since  ©i*  is  nonsingular,  we  have 

U[N,H]=0.  (A.  17) 

In  Section  4.1.1  we  showed  that  the  (i,  j)th  element  of  the  skew-symmetric  matrix 
[N,  H ]  is  hij  ( rii  —  rij).  Equation  A. 17  tells  us  that  the  (1,  2)th  and  (3, 4)th  elements 
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of  [TV,  H]  are  zero.  Since  the  elements  of  TV  are  distinct,  we  have 


h\2  —  h2 1  —  0 

(A- 

h,34  =  h.43  =  0 

.  Hence,  the  gradient  flow  equation  for  the  1-level  ROT  is  guaranteed  to  converge 
to  a  ©i*  such  that  the  resulting  approximately  diagonalized  matrix  H  has  h\2  = 

Ti-34  =  o. 


A.3.2  2-Level  ROT 

For  the  2-level  ROT  we  have 

0  =  0^202,  (A. 19) 

where  ©i  and  ©2  are  of  the  form  given  in  Equation  A.  2  and  P2  is  given  by  Equa¬ 
tion  A. 3.  Following  Section  4.1.2,  the  gradient  flow  equations  to  search  for  the 
ROT  variables  which  most  nearly  diagonalize  the  symmetric  matrix  H0  are 

©i  =  -©in  [P2©2TV0pf,  ©ftfoBi]  ,  0i(O)  =  010 

02  =  -02n  [TV,  ©^©rPo01P202]  ,  ©2(0)  =  ©20-  (A. 20) 

Let  (©i* ,  ©2*)  be  an  equilibrium  point  for  this  system  of  equations  and  let 

H  =  ©£P2t  ®IH0QuP2Q2*  (A. 21) 

be  the  resulting  approximately  diagonalized  matrix.  Looking  at  the  equation  for 
©2  and  following  the  calculations  from  the  previous  section,  we  see  that  hi2  = 
h2\  =  h34  =  /i43  =  0.  Looking  at  the  Equation  for  ©4,  we  see  that 

n  [P202*TV0£P2t,  ©[PoBi]  =  0.  (A. 22) 


138 


Using  some  algebraic  manipulation  we  have 


[p2©2*iv©^p2T,  e^H0eu] 

=  [P202*lVe£P2T,  P202*©£P2T ©uPo©i*P2©2*©^Pj] 

=  P202*  [N,  H]  &lPj.  (A. 23) 


We  can  write  02*  as  a  matrix  using  sines  and  cosines, 

cos($2i)  —sin(' #2i)  0  0 

sin(t?2i)  cos(i?2i)  0  0 

0  0  cos(t9  22)  -sin(i9  22) 

0  0  sm(i?  22)  cos($  22) 

.  To  simplify  upcoming  long  expressions,  we  will  use  and  Qj  to  denote  sin('dij) 
and  cos(djj),  respectively. 

Now,  the  equilibrium  condition  means  that  the  (1,  2)th  and  (3, 4)th  elements  of 
the  matrix  P202*  [N,  H]  0^P2T  must  be  equal  to  zero.  Substituting  Equation  A. 24 
into  this  expression  and  evaluating  using  the  Maple  kernel  of  MATLAB  yields  the 
following  two  equations: 

^13s22C2l(u.3  —  n  1)  +  h14C22C2l(n4  —  rii) 

+h23s22s21(n2  -  n3)  +  h2Ac22s21(n2  -  n4)  =  0 

( A.  Z10  ) 

hi3c22s2i(ri3  —  rii)  +  /i14S22s2i(n4  —  u4) 

+^23c22c21  {p>2  —  n3 )  +  h2AS22C2i(n2  —  U4)  =  0. 

Hence,  the  gradient  flow  equations  for  the  2-level  ROT  are  guaranteed  to  con¬ 
verge  to  an  equilibrium  point  (©1*,02*)  such  that  the  resulting  approximately 
diagonalized  matrix  H  has 


1.  h\2  —  h^4  —  0 
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2.  The  vector  [/ii3/ii4/i23/i24]T  is  in  the  null  space  of  the  matrix 

S22C2i(n3  -  ni)  c22c21(nA  -  m)  s22s21(n2  -  n3)  c22s2i(n2  -  n4) 
c22s2i(n3  —  ni)  s22s2i(n4  —  n4)  c22c2i(n2  —  n3)  s22c2i(n2  —  n4) 

A. 4  3— Level  ROT 


Here  we  consider  the  3-level  ROT  given  by  Equation  A.l.  Following  Section  4.1.2, 
the  gradient  flow  equations  to  search  for  the  ROT  variables  which  most  nearly 
diagonalize  the  symmetric  matrix  //()  are 

©i  =  -©in  [P202P303iV0jP3r0^P2T,  ©fifoBi]  ,  0i(O)  =  010  (A. 26) 

02  =  -02n  [P303iV0[P3T,  0^0^00^202]  ,  02(O)  =  020  (A. 27) 

03  =  -03n  [TV,  0^P3T0^P2T0fPo0in202P303]  ,  ©3(0)  =  030-  (A. 28) 

Let  (0i*,  02*,  ©3*)  be  an  equilibrium  point  for  this  system  of  equations  and  let 

H  =  ©£P3T©^©f*tfo©l*P2©2*P3©3*  (A. 29) 

be  the  resulting  approximately  diagonalized  matrix.  Looking  at  the  equation  for 
03  and  following  the  calculations  from  the  Section  A. 3.1,  we  see  that  /i12  =  h2 1  = 
h34  =  h43  =  0.  Looking  at  the  equation  for  02  and  following  t  the  calculations  of 
Section  A. 3. 2  we  get  the  equations 

hi3s32c31(n3  -  n4)  +  h14c32c3i(n4  -  n4) 

+h23S32S3i(^2  -  n3)  +  h24c32s3i(n2  -  n4)  =  0 

^i3c32-s3i(n3  -  m)  +  h14s32s31(n4  -  n4) 

+h23c32C3i(n2  -  n3)  +  h2As32c3i{n2  -  n4)  =  0, 
where  s3j  and  c3;  are  the  sines  and  cosines  that  result  from  writing  03*  as  a  matrix 

of  sines  and  cosines  (see  Equation  A. 24). 


(A. 31) 
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(A. 32) 


Now,  looking  at  the  Equation  for  ©i,  we  see  that 

n  [p3©3*p2©2*tv©£p2t  ©£p3t,  elH0ou]  =  o. 

Using  some  algebraic  manipulation  we  see  that 

[P202*iV0^P2T,  0f*po0i*]  =  P202*P303*  [N,  H]  ©r*P3T©£P2T  (A. 33) 

The  equilibrium  condition  means  that  the  (1,  2)th  and  (3, 4)th  elements  of  the  ma¬ 
trix  P202*P303*  [N,  H]  QIPTQIP?  musf  be  equal  to  zero.  Using  the  sine/cosine 
matrix  representations  of  02*  and  03*  along  with  the  Maple  kernel  of  MATLAB 
yields  the  equations 

^i3(-s2is32s3ic22  +  c2ic32c3is22)(n3  —  ni) 

+/ii4(s2iC32s3iC22  —  c2is32c3is22)(n4  —  rii) 

+h23(s2is32c3ic22  —  c2ic32s3is22)(n3  —  n2) 

+/i24(s2ic32c3ic22  +  c2is32s3is22)(n4  —  n2)  =  0 


(A. 34) 


(A. 35) 


hi3(c2is32s3is22  +  s2ic32c3ic22)(ni  —  n3) 

+hi4(c2ic32s3is22  —  S2i-s32c3ic22)(ni  —  n4) 

+  ^23(c2lS32C3lS22  —  S2lC32S3iC22)  (n2  —  U3) 

+h24(c2ic32c3is22  +  s2i-s32s3ic22)(n2  —  n4)  =  0. 

Since  hi3,  hi4,  /i23,  and  /i24  all  enter  linearly  into  Equations  A. 30,  A. 31,  A. 34, 

and  A. 35,  we  can  rewrite  these  equations  as 


(A. 36) 


where  T  is  a  4x4  matrix.  When  the  matrix  T  is  full  rank,  then  we  must  have 
hi3  =  hi4  =  h23  =  h2 4  =  0  at  the  equilibrium  points  of  the  ROT  gradient  flow 


CO 

0^ 

0 

hi4 

0 

^23 

0 

h2  4 

0 
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equations.  From  the  03  equation  we  already  have  that  h  1 2  =  /i34  =  0,  So  when  T 
is  full  rank,  the  matrix  H  must  be  diagonal. 

Using  the  Maple  kernel  of  MATLAB,  we  see  that  the  determinant  of  T  is 

det(r)  =  (n4  -  rq)(n2  -  n4)(n3  -  rq)(n3  -  n2)(c222s22l  -  s222c221).  (A.37) 

Hence,  the  determinant  is  nonzero  and  I’  is  full  rank  unless 

cos2('t?22)sin2(t?2i)  =  sin2  ($22 )  cos2  ($21 )  •  (A.38) 

This  failure  in  rank  leads  to  the  possibility  that  the  matrix  H  has  nonzero  off- 
diagonal  elements. 

We  can  take  some  comfort  in  the  fact  that  T  is  nonsingular  “almost  every¬ 
where”,  i.e.  the  set  on  which  T  is  singular  is  a  thin  subset  of  the  configuration 
space  of  the  ROT  variables.  It  is  tempting  to  state  that,  for  most  cases,  the 
gradient  flow  equations  for  the  3-level  ROT  to  converge  to  an  equilibrium  point 
(©1*,  ©2*,  ©3*)  such  that  the  resulting  H  is  diagonal.  Based  on  our  numerical  ex¬ 
periments,  this  fact  seems  to  be  true.  However  the  theoretical  argument  outlined 
above  is  incomplete.  It  is  entirely  possible  that  the  thin  set  on  which  T  is  singular  is 
an  attractor  for  the  system  of  ROT  equations.  In  order  to  complete  our  argument, 
it  is  necessary  to  rule  out  the  possibility  that  the  ROT  equations  converge  to  an 
equilibrium  point  which  lies  on  this  thin  set.  This  remains  as  an  open  problem. 
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Appendix  B 


Flexible  Beam  Model 


This  appendix  develops  a  model  for  a  flexible  cantilever  beam.  The  beam  is  driven 
by  a  number  of  surface  mounted  PZT  actuators  and  the  beam  has  a  number 
of  outputs  given  by  surface  mounted  PZT  sensors.  We  assume  that  the  beam 
undergoes  pure  bending  and  satisfies  the  Euler-Bernoulli  condition  (linear  strain 
distribution).  From  these  assumptions,  we  use  first  principles  to  derive  a  PDE 
model  to  describe  the  dynamics  of  the  beam.  The  PDE  is  then  discretized  using 
the  finite  difference  method  to  give  an  ODE  model.  Simulation  results  using  this 
ODE  model  are  then  presented. 

B.l  Simple  Beam  Model 

Here  we  derive  the  equations  of  motion  for  transverse  vibrations  in  a  cantilever 
beam.  The  beam  has  length  L,  width  b,  and  thickness  4-  We  assume  that  L 
b  tc.  At  rest,  the  beam  is  oriented  as  depicted  in  Figure  B.l.  The  clamped  end 
of  the  beam  is  at  x  =  0  and  the  free  end  is  at  x  =  L. 

It  is  necessary  to  introduce  some  notation  with  which  we  can  describe  the  beam 
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Figure  B.l:  Rest  configuration  of  beam. 

when  it  is  deformed  in  some  way.  One  common  method  of  doing  this  is  to  describe 
the  displacement  of  each  point  in  the  beam  using  rest  position  as  a  reference  point. 
Since  there  are  a  continuum  of  points  within  the  beam,  these  descriptions  take  the 
form  of  continuous  functions  in  x,  y,  and  z.  In  order  to  describe  how  the  beam 
evolves  over  time,  these  functions  must  also  be  time  dependent.  Let  x,  y,  and  z 
be  the  coordinates  of  a  point  P  in  the  beam  at  rest.  Define  the  x  component  of 
the  displacement  of  the  point  P  at  time  t  by  u(x,y,z,t).  Likewise,  let  the  y  and 
z  components  be  denoted  by  v(x,  y,  z,  t )  and  w(x,  y,  z,  t ),  respectively. 

For  the  current  case,  only  transverse  vibrations  are  considered.  We  assume 
that  there  is  no  displacement  in  the  y  direction  which  means  that  v  is  identically 
zero.  We  assume  that  the  beam  does  not  undergo  extension  which  means  that 
u  is  identically  zero.  We  also  assume  that  there  is  no  torsion  in  the  beam,  so  w 
is  not  dependent  on  y.  Since  the  beam  is  thin,  we  assume  that  the  displacement 
function  does  not  vary  across  the  thickness  of  the  beam.  As  a  result  of  these 
assumptions,  the  configuration  of  the  beam  can  be  completely  described  using  the 
function  w  =  w(x,t),  as  shown  in  Figure  B.2. 

At  time  t,  we  define  the  neutral  surface  to  be  {z\z  =  w(x,t),x  G  [0,L]}.  The 
strain  in  the  neutral  surface  is  zero  as  a  result  of  the  “no  extension”  assumption. 
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Figure  B.2:  Beam  displacement. 


This  is  shown  in  Figure  B.2.  Away  from  the  neutral  surface,  the  strain  is  assumed 
to  be  axial;  shear  strain  is  neglected. 

We  also  assume  that  the  beam  is  an  Euler-Bernoulli  beam.  This  assumption 
means  that  each  cross  section  perpendicular  to  the  x  axis  for  the  resting  beam 
remains  perpendicular  to  the  neutral  surface  for  all  time.  This  assumption  can  be 
used  to  find  a  formula  for  the  strain  at  any  point  in  the  beam.  Figure  B.3  depicts  a 
differential  element  of  the  beam  undergoing  some  deformation.  The  cross  sections 
at  x  and  x  I  ■  dx  are  perpendicular  to  the  neutral  surface,  which  is  shown  here  as 
an  arc  connecting  the  two  points.  The  angle  between  the  two  cross  sections,  9,  is 
shown  to  be  equal  to  the  negative  curvature  of  the  beam, 

d2w 


6  = 


dx 2 


(x)dx. 


(B-l) 


The  radius  of  curvature,  R ,  can  be  found  using  the  relationship 

dx  9  1  d2w  ,  . 

— ^  =  —  = - 7rxx(x)dx, 

2nR  2-7T  2n  dx 2 

which  leads  to  the  following  expression  for  R: 

'  d2w , 


(B.2) 


R  = 


dx2 


[x, 


-i 


(B.3) 
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Figure  B.3:  Kinematic  diagram  of  differential  beam  element  used  to  determine 
strain. 


To  calculate  the  strain  at  a  point  on  the  cross  section,  it  is  necessary  to  look  at 
the  change  in  length  of  the  “fiber”  passing  through  that  point.  A  fiber  is  defined 
to  be  a  set  of  points  in  the  resting  beam  which  form  a  line  parallel  to  the  x  axis. 
When  the  beam  is  deformed,  the  fibers  are  also  deformed.  In  Figure  B.3,  the  arc 
labeled  “fiber  section”  represents  the  piece  of  a  deformed  fiber  which  is  contained 
in  the  differential  element.  The  fiber  is  displaced  from  the  neutral  surface  by  £. 
The  length  of  the  arc  is  denoted  by  a{Q.  In  the  undeformed  beam,  the  length  of 


the  fiber  section  is  dx ,  which  is  the  same  as  the  length  of  the  differential  element. 
In  the  deformed  beam,  we  can  obtain  an  expression  for  a(()  using  the  relationship 


dh  =±  =  -±^ix)dx 

(R  +  0  2x  dx 2 


(B.4) 


which  yields 


a(C)  =  dx-  ( 


d2w 

dx2 


(x)dx. 


(B.5) 


The  strain  in  the  fiber  section,  which  is  defined  to  be  the  change  in  length  divided 
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by  the  undeformed  length,  is  given  by 


6(C) 


a(C)  —  dx 
dx 
d2w 


Therefore,  at  time  t  the  strain  at  a  point  P  in  the  perpendicular  cross  section  at 
x  is  normal  to  the  cross  section  with  magnitude 


e(C^)  =  -C^Jr(:M),  (B.6) 

where  Q  is  the  displacement  between  P  and  the  neutral  surface. 

To  get  the  relationship  between  stress  and  strain,  we  use  a  version  of  Hooke’s 
law  modified  to  account  for  viscous  damping  within  the  beam  [4] .  This  law  yields 
the  stress  on  the  cross  section  as 


MO  =  M  + 

=  (B.T) 

where  Eb  and  E£  are,  respectively,  Young’s  modulus  and  the  viscous  damping 
constant  for  the  beam  material.  Here  the  subscript  b  emphasizes  that  ab,  Eb,  and 
El  correspond  to  the  beam.  The  reason  for  this  emphasis  will  become  apparent  in 
Section  B.2.  The  force  acting  on  a  differential  area  dA  of  the  beam  cross  section 
is  then  crb(C)dA.  This  is  shown  in  Figure  B.4.  Here  we  have  chosen  the  convention 
that  a  normal  force  directed  away  from  the  cross  section  is  positive. 

The  bending  moment  acting  on  the  cross  section  can  now  be  computed. 


Mb(x,t ) 


-  /  ^(C)C  dA 

J  A 

PI  2 


-4/2 


ab(()b(d( 
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Figure  B.4:  Strain  in  beam  cross  section. 


ft/2  (  de 

Lb,AEbe+E^]KdC 


ft/2 


=  (Et 


-4/2 

d2w 


02  w  d  02  w 

^c^+^(cS)mc 


dx 2 


+  e : 


d3w  , 
dtdx 2 ' 


dt^dx2 

ft/2 


'-4/2 


_  r  <92te  <93w 

Ebh-TTT  +  EbIb- 


dx 2 


<9tdx2’ 


(B-8) 


where  Ib  is  the  moment  of  inertia  of  the  cross  section.  In  this  case,  //,  =  bt\jYl. 

Now  we  consider  the  forces  and  moments  acting  on  an  interior  differential 
element  of  the  beam  as  shown  in  Figure  B.5.  Balancing  the  forces,  we  get 


n  ,n  dS  d2w 

S~(S+  ~  Pbtbdx~w 


=  0 


pbtb 


d 2 


w 


dt2 


dS_ 

dx 


(B-9) 


Ignoring  the  rotational  inertia  of  the  differential  element,  the  moment  balance 
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There  is  also  no  shear  force  at  the  free  end,  which  implies 


EAdE(L't)*E’lhMd7*(L’t)  =  0  Vi>0'  (B'16) 

Given  initial  conditions  w(x,  0)  and  j^r(x,  0)  for  x  G  [0,  L\,  Equations  B.12  through 
B.16  give  a  complete  description  of  the  deformation  of  the  beam  for  all  t  >  0.  1 


B.2  Sandwiched  Beam  Model 

In  this  section,  we  derive  the  equations  of  motion  for  transverse  vibrations  in  a 
flexible  beam  which  is  sandwiched  between  two  PZT  actuators.  The  actuators 
are  mounted  with  opposite  polarity  so  that  applying  the  same  voltage  to  both 
actuators  induces  pure  bending  (no  extension)  in  the  beam. 

A  voltage  applied  across  the  terminals  of  a  PZT  actuator  induces  a  stress 
tensor  field  within  the  crystal,  which  in  turn  induces  a  strain  in  the  crystal.  In  this 
model,  we  consider  only  the  component  of  stress  in  the  axial  direction;  all  other 
components  are  ignored.  An  applied  voltage  Vc  will  create  an  internal  stress  of 

= -d31Ec^  +  Ecec,  (B.19) 

ic 

in  the  top  crystal  and 

abott°m  =  d3iEVc  +  (B.20) 

lc 

1If  the  damping  is  small  compared  to  the  stiffness  (Ef,  E£),  boundary  conditions  B.15  and 
B.16  can  be  approximated  by  the  simplified  boundary  conditions 

=  o  (B.17) 

and 

5f(C<)=  0,  (B.  IS) 

respectively,  for  allf  >  0. 


150 


Figure  B.6:  Strain  in  sandwich  cross  section. 

in  the  bottom  crystal  [13].  In  these  equations  d31  is  the  piezoelectric  strain  coeffi¬ 
cient,  Ec  is  Young’s  modulus  for  the  PZT  crystal,  tc  is  the  thickness  of  the  PZT, 
and  ec  is  the  axial  strain  in  the  PZT.  We  neglect  the  viscous  damping  term  in  the 
crystal  for  simplicity,  we  will  keep  the  viscous  damping  term  in  the  beam  material. 
The  parameter  d3 1  is  assumed  to  be  negative  so  that  a  positive  voltage  causes  the 
top  crystal  to  contract  while  the  bottom  expands,  creating  a  positive  curvature. 

We  assume  that  the  Euler-Bernoulli  bending  assumption  holds  in  the  crystal  as 
well  as  the  beam,  so  the  strain  in  the  crystal  is  a  linear  function  of  beam  curvature 
and  distance  from  the  neutral  surface, 

d^w 

^(0  =  -CW.  (B.21) 

Now  we  can  compute  the  bending  moment  acting  on  the  sandwich  cross  section. 
The  forces  used  to  compute  the  moment  are  shown  in  Figure  B.6. 

Ms(x,t ) 
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(B.22) 


where  Mb  is  given  by  Equation  B.8  and  /,  is  the  moment  of  inertia  of  the  crystal 
cross  sections  about  the  neutral  axis.  In  this  case,  Ic  =  {b{2tc  E  b>)3  —  bt\)/ 12.  2 
As  in  the  case  of  the  simple  beam,  the  forces  and  moments  balance  to  give 

d2w  d2Mq 


(2 Pete  E  Pb^b)  b~ 


dt 2 


dx 2 


(B.23) 


Substituting  the  expression  for  Ms  from  Equation  B.22  into  Equation  B.23  yields 

the  equation  of  motion  for  the  interior  of  the  sandwich. 

d2w  d5w 

(2 Poto  +  pA) bw  =  -  (ECIC  +  EtIb)  -  Elhg^.  (B.24) 

Note  that  the  forcing  term  in  Ms  due  to  the  PZT  voltage  is  constant  in  x  and 
disappears  when  differentiated.  This  term  only  affects  the  boundary  elements  of 
the  sandwich.  Specifically,  at  x  =  L  there  is  no  externally  applied  moment,  which 
means 

d2w  d3w 

bd31Ec(tb  E  tc)Vc  +  (ECIC  E  EhIb)  E  E*bIb-^(L,t)  =  0  (B.25) 


2Here,  we  have  assumed  that  the  crystal  thickness  tc  constant  and  that  the  width  of  the 
crystal  is  everywhere  equal  to  the  width  of  the  beam.  In  some  applications,  it  may  be  desirable 
to  allow  one  or  both  of  these  parameters  to  vary,  causing  Is  and  the  number  which  multiplies  Vc 
to  be  functions  of  x. 
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for  all  time  t  >  0.  The  shear  force  at  x  =  L  is  also  zero,  which  implies 


(ECIC  +  EbIb)  t)  +  Ethg^(L,  t)  =  0  (B.26) 

for  all  time  t  >  0.  At  the  clamped  end,  the  position  and  slope  of  the  beam  are 
both  zero,  yielding 


w(0,t)  =  0  (B.27) 

£f«M)  =  0  (B.28) 

(B.29) 

for  all  time  t  >  0.  Given  initial  conditions  w(x,  0)  and  ^-(x,0)  for  x  G  [0,  L], 
Equation  B.24  and  Equations  B.25  through  B.28  give  a  complete  description  of 
the  deformation  of  the  beam  for  all  time  t  >  0.  3 


B.3  Final  Beam  Model 


In  this  section,  we  derive  the  equations  of  motion  for  a  flexible  cantilever  beam 
excited  by  pairs  of  PZT  crystals  mounted  at  various  places  on  the  beam.  We 
assume  that  all  of  the  crystals  have  the  same  thickness,  tc.  Additionally,  the 


3If  the  damping  coefficient  of  the  beam  is  small  compared  to  the  stiffnesses  of  the  beam  and 
the  crystal  (i.e.  is  much  smaller  than  Ef,  and  Ec),  then  the  boundary  conditions  given  by 
Equations  B.25  and  B.26  can  be  approximated  by 


d2u>  bd31Ec(tb  +  tc)Vc 

dx*[  ’  j  “  ECIC  +  EbIb 


(B.30) 


and 


respectively,  for  allt  >  0. 


d3w 

dx3 


=  0, 


(B.31) 
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x=X1  x=ji|  x=X2  x=jj.2  x=A,3 

x=0 


x=L 


Figure  B.7:  Example  configuration  of  final  beam. 

crystals  all  have  the  same  width,  b ,  which  is  also  the  width  of  the  beam.  An 
example  of  such  a  beam  is  depicted  in  Figure  B.7. 

First,  we  introduce  some  notation  for  describing  the  locations  of  the  PZT  crys¬ 
tals.  Let  k  be  the  number  of  PZT  pairs.  We  will  use  A*  to  denote  the  x  coordinate 
of  the  left  edge  of  the  zth  PZT  pair.  Let  /q  denote  the  x  coordinate  of  the  right 
edge  of  the  ith  PZT  pair.  The  PZTs  cannot  overlap  and  they  cannot  extend  be¬ 
yond  the  dimensions  of  the  beam.  We  assume  that  the  first  pair  is  the  pair  closest 
to  the  left  end  of  the  beam.  We  also  define  fj, o  =  0  and  A^+i  =  L.  Then  we  have 

0  =  ^  Ai  <  /xi  <  A2  <  <  •  •  •  <  Afc  <  //fc  <  Afc+i  =  L.  (B.32) 

We  will  call  the  sections  of  the  beam  with  no  PZT  “simple”  sections  and  we  will 
call  the  sections  with  PZTs  “sandwiched”  sections. 

Within  the  simple  sections,  the  motion  of  beam  is  described  by  Equation  B.12. 
Likewise,  Equation  B.24  describes  the  motion  of  the  interiors  of  the  sandwiched 
sections.  In  order  to  obtain  a  complete  model  for  the  beam,  it  is  necessary  to 
piece  these  sections  together  with  the  appropriate  boundary  conditions.  The  left 
end  of  the  beam  is  clamped,  so  the  boundary  conditions  at  x  =  0  are  given  by 
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Equations  B.13  and  B.14.  At  x  =  L,  the  beam  is  free  so  the  boundary  conditions 
are  given  by  Equations  B.15  and  B.16  if  the  rightmost  section  of  the  beam  is 
simple  or  Equations  B.25  and  B.26  if  the  rightmost  section  is  sandwiched.  We 
have  yet  to  determine  the  boundary  conditions  at  the  points  where  simple  section 
meets  sandwiched  section  (i.e.  x  G  {Ai,  fi\, . . . ,  A Hk})-  Let  V*  denote  the  voltage 
applied  to  the  zth  PZT  pair. 

First  we  need  to  introduce  some  more  notation.  Consider  the  function  /  :  R  — > 
R.  We  define  the  /(x q  )  to  be  the  limit  of  /(x)  as  x  approaches  Xo  from  the  left, 

i.e. 

f(xo)  =  i™  fix).  (B.33) 

X<XQ 

Similarly,  we  define  /(xq)  to  be  the  limit  of  /(x)  as  x  approaches  x0  from  the 
right. 

We  assume  that  w  and  ^  are  continuous  functions  in  x,  which  means  that  for 
all  t  >  0 


^(A  i,t)  = 

(B.34) 

dw 

(B.35) 

at  the  left  sandwich-beam  boundaries  and 

(B.36) 

dw , 

= 

dw .  ,  . 

(B.37) 

at  the  right  boundaries.  We  assume  the  bending  moment  is  continuous  in  the 
beam,  so 

M(K)  =  M(\ t).  (B.38) 


Equations  B.8  and  B.22  are  substituted  into  Equation  B.38  to  yield 


EbIb—(\-,t)  =  bd31Ec(tb  +  tc)v;  +  (EbIb  +  ECIC)  ^-(A+t). 


dx2 


(B.39) 
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This  gives  us  one  set  of  boundary  conditions  for  the  left  edges  of  the  PZT.  The 
same  process  computed  at  p,,;  gives  a  set  of  boundary  conditions  for  the  right  edges 
of  the  PZTs, 


d2w  d2w 

bd3iEc(tb  +  tc)V*  +  ( EbIb  +  ECIC )  1 1)  —  >  ^)-  (B.40) 


Similarly,  we  assume  the  shear  force  is  continuous  in  the  beam.  Setting  S(x  )  = 
S(x+ )  at  the  boundaries  yields  additional  boundary  conditions 


EbIb^(K,t)  =  (E„h  +  ECIC)  |^(A  +  t) 


(B-41) 


(E„h  +  ECIC)  t)  =  EbIb  ^(j4,  t). 


(B.42) 


B.3.1  Final  Model  Statement 


We  can  use  the  above  information  to  piece  together  a  model  for  the  entire  beam. 
It  is  stated  here: 


simple  sections: 

For  x  e  (Hi-i,  \),  i  =  {1, 2,  3, . . . ,  k  +  1}, 


,  d2w  d4w  „  d5w 

pbtbW  =  ~  bhd %  “  bhdtd. 


(B.43) 


sandwiched  sections: 

For  x  e  (A t,  fJ,i),i  =  {1,  2,  3, ... ,  k}, 

d2w  d4w  d5w 

<**•  +  pA)  b—  =  -  (EX  +  ^  (B.44) 

boundary  conditions  (clamped  end): 


w(0,t)  =  0 


(B.45) 


-(0  ,t)  =  0 


(B.46) 


for  all  time  t  >  0. 


boundary  conditions  (free  end): 

If  the  rightmost  section  of  the  beam  is  simple,  then  for  all  t  >  0 


d^w  d^w 

=  o 

d^w  d^w 

Ebh—(L,t)  +  Elh—(L,t)  =  o. 


(B.47) 

(B.48) 


If  the  rightmost  section  of  the  beam  is  sandwiched,  then  for  all  t  >  0 


bd31Ec(tb  +  tc)Vck{t )  +  (ECIC  +  EbIb)  gf  ( L,t )  +  E*bIb^(L,t) 


(B.49) 


(EJC  +  EA)  §^(L,  t )  +  Eth-^iL,  0  =  0. 

boundary  conditions  (left  edges  of  crystals): 


(B.50) 


=  w(xi^) 


dw  , ,  _  , 


■(V.  0 


(B.51) 

(B.52) 


M3i^c(4  +  *C)V£  +  (^4  +  ECIC)  —{X f,  t )  (B.53) 


(EbIb  +  EcIc)^(Xt,t) 


(B.54) 


for  i  =  {1, 2, . . . ,  /c},  Vt  >  0. 

boundary  conditions  (right  edges  of  crystals): 


dw  , 
d^w 

bd:nEc(tb  +  tc)V*  +  (EbIb  +  ECIC )  2  (/j,i  ,  t) 

crm 

(EbIbEEcI^—(p-A) 


(B.55) 


(B.56) 

E,lll—<i4.l>  (B.57) 
(B.58) 
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for  i  =  {1,  2,  3 , ,k},  Vf  >  0  if  the  rightmost  section  of  the  beam  is  simple,  for 
i  =  {1,2,3,  — 1},  V*  >  0  if  the  rightmost  section  of  the  beam  is  sandwiched. 

Given  initial  conditions  w(x,  0)  and  ^f(x,  0)  for  x  G  [0,  L\,  and  inputs  V*(t),  i  = 
{1,2,...,  k},  Vf  >  0,  these  equations  give  a  complete  description  of  the  deformation 
of  the  beam  for  all  time  t  >  0. 


B.4  Finite  Difference  Model 


Here  we  apply  the  finite  difference  method  to  the  partial  differential  equations 
(PDEs)  and  boundary  conditions  listed  in  Section  B.3.1  to  obtain  a  set  of  ordinary 
differential  equations  (ODEs)  which  approximate  the  PDE  system. 

The  first  step  in  applying  the  finite  difference  method  is  to  discretize  the  spatial 
domain.  For  the  beam,  we  consider  N  evenly  spaced  points  between  x  =  0  and 
x  =  L.  We  will  include  the  point  at  x  —  L.  We  do  not  include  the  point  as  x  =  0 
(it  is  not  necessary  to  include  this  point  since  the  boundary  conditions  tell  us  is 
it  always  zero.)  We  label  the  leftmost  point  X\  and  the  rightmost  point  xjy.  This 
allows  us  to  write  the  zth  point  as  Xi  =  iA,  where  A  =  L/N  is  the  spacing  between 
the  points.  We  denote  the  displacement  of  the  beam  at  the  point  xt  and  time  t  as 

Wi(t )  =  w(xi,t).  (B.59) 


The  next  step  is  to  approximate  the  spatial  derivatives  of  w  using  the  w,s.  First 
some  notation: 

r)^ii l-  a  rfim  l 

(B.60) 


dxk  dxk 

The  approximations  for  the  first  four  spatial  derivatives  are  then 


dwi  Wi  -  Wi- 1  dwi  wi+1  -  Wi 

dx  ~  A  °r  dx  ~  A 


(B.61) 
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(B.62) 


d 2 


Wj 


Wi-i  -  2 Wi  +  wi+1 


dx2 


d3 


Wj 


dx3 


-Wi- 2  +  3wi-i  -  3^i  +  wi+i 

A3 


or 


A2 

d3Wi 

dx3 


-w^ i  +  3 Wi  -  3wi+i  +  wi+2 


A3 


and 


d4Wi  Wi- 2  -  4u>j_i  +  6uy  -  Awi+1  +  wi+2 


(B.63) 


(B.64) 


dxA  A4 

Of  course,  Wi  is  only  defined  for  i  —  { 1,2,...,  N},  so  these  approximations  are  only 

valid  at  the  “interior  points”,  which  are  loosely  defined  here  as  points  far  enough 
away  from  the  boundaries  that  the  above  spatial  derivative  approximations  are 
valid. 

We  will  obtain  a  system  of  ODEs  approximating  the  PDE  system  by  substitut¬ 
ing  Equations  B.61,  B.62,  B.63,  and  B.64  into  Equations  B.43  and  B.44  for  each 
Wi .  First,  we  introduce  some  more  notation: 

A  d 


Wi  = 


Wj  = 


Wj 


rfe 

■p* 

1  b 

Tr 


—Wi 

dt 

dt2 
Eblb 
Pbbtb 

E*bIb 

pbtb 

Eblb  +  ECIC 


(B.65) 

(B.66) 

(B.67) 

(B.68) 

(B.69) 


b(pbtb  +  2  pctc) 

Now  we  can  use  this  notation  to  rewrite  the  equations  of  motion  for  the  interior 
of  a  simple  beam  section: 

_  d5w 


d2w  d4w 

dt2  b  dxA 

Subsituting  from  Equation  B.64  yields 

"i  i  "p*  d 


’  dt  dxA ' 


(B.70) 


Wi=  - 


A4 


(Wi- 2  -  4iUj_i  +  6 Wi  -  4wi+i  +  wi+2) 


(B.71) 
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Again,  this  expresseion  is  valid  only  at  the  “interior  points”  in  the  simple  section, 
i.e.  points  Xi  such  that  x*_2,  Xj_i,  x*,  Xj+i,  and  Xi+ 2  are  all  contained  in  the  simple 
section.  Hence,  Equation  B.71  provides  a  second  order  ODE  for  each  interior  point 
in  the  simple  section.  Likewise,  within  a  sandwiched  section  we  get 

(vc  +  r*— \ 

Wi  =  -  I  c  ^4cdt  |  (wi- 2  -  Awi-x  +  Qwi  -  4wi+1  +  wi+2)  ■  (B.72) 

Again,  this  expression  is  only  valid  for  interior  points  of  the  sandwiched  section. 

So  we  have  now  derived  second  order  ODEs  to  describe  the  motion  of  each 
interior  point  on  the  beam.  Now  we  need  to  find  ODEs  to  describe  the  remaining 
“boundary  points”.  This  can  be  done  using  the  boundary  conditions  from  the 
PDE  model. 

The  problem  that  arises  at  the  boundary  points  is  that  the  spatial  derivative 
approximations  require  the  values  of  displacements  at  points  that  are  not  defined 
for  the  original  PDE.  To  solve  this  problem,  we  define  “ghost  points”.  The  dis¬ 
placement  at  these  ghost  points  can  be  found  using  the  boundary  conditions. 

For  example,  at  the  clamped  end  of  the  beam,  we  define  ghost  points  .x0  and 
X-\  which  lie  A  and  2 A  to  the  left  of  X\ ,  respectively.  We  denote  the  displacement 
at  Xi  by  Wi.  The  boundary  condition 

w(0,t)  =  0  Vt>0  (B.73) 


implies  that  w0  =  0  Vf  >  0.  The  boundary  condition 

dw 

V(>0 

implies 

w0  -w_  1  _ 

A 

=>•  w-i  =  0  Vf  >  0. 


(B.74) 


(B.75) 


160 


Using  these  ghost  points,  we  can  now  write  a  valid  expressions  for  w\  and  ib2 
(assuming  the  leftmost  section  of  the  beam  is  simple): 

^  +  Tb-dt 


Wi  =  - 


(ic_i  —  4  w0  +  6?a!i  —  4w2  +  w3) 


(6u!i  —  4u)2  +  w3)  (B.76) 


and 


w2  =  - 


(wq  —  4wi  +  6w2  —  4te3  +  W4) 


(— 4wi  +  6«;2  —  4w3  +  iu4) .  (B.77) 

At  the  free  end  of  the  beam,  we  introduce  ghost  points  x^+i  and  x^+2.  As¬ 
suming  the  rightmost  section  of  the  beam  is  simple,  the  zero  moment  condition 


is 


^(L,t)  +  £.4^_(L,t)=0  V(>0. 

Subsituting  the  appropriate  spatial  derivative  approximations  yields 

wn-i  —  2wn  +  ^wn+i  \ 


(B.78) 


Ebh 


(  %- 1  —  2wn  +  Wn+ 1  \ 
1  ^  ) 


E*bh 


A2 


=  0,  (B.79) 


which  simplifies  to  the  ODE 


d EbIb  /  ~  \ 

37 WN+1  =  -7 777-  (%-l  -  2wat  +  %+l)  -  Wjv-1  +  2Wiv, 
at  ^bE 

which  must  hold  for  all  t  >  0.  The  shear  condition  at  the  free  end  is 

^(£-*)  +  £I^(£-*)  =  0  Vi>0' 

Making  the  spatial  derivative  approximations  yields 

1  +  3^  -  3tCAT+i  +  wat+2\ 


(B.80) 


(B.81) 


Ebh 


+E*bh 


A3  J 

—Wn-1  +  3wn  —  3|wat+i  +  ^Wn+2 
A3 


=  0. 


(B.82) 
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Simplifying  and  substituting  for  ^jWn+i  from  Equation  B.80  yields  the  ODE 


d  _  Eblb  _  . 

3Tu’Ar+2  =  ~  pz  T  (2u)7V— 1  —  3wat  +  WjV+2)  —  2wjv-1  +  3u!jv,  (B.83) 

at  Eh  If, 

which  hold  for  all  t  >  0. 

Using  these  ghost  points,  we  can  now  write  a  valid  expressions  for  Wn-i  and 
Wn  (assuming  the  rightmost  section  of  the  beam  is  simple): 

■>  _|_  p*_d_ 

WN~  1  =  -  1  - A  Abdt  |  (wN-3  -  4u)jv— 2  +  6wN-i  -  4wn  +  wN+i)  (B.84) 


and 


wN  = 


A4 


+  Tbdt 

A4 


(wN- 2  -  4wat_i  +  6 wN  -  4wn+1  +  wN+2) ,  (B.85) 


where  the  quantities  vjn+i  and  Wn+2  are  given  by  the  solutions  of  the  ODEs  in 
Equations  B.80  and  B.83  respectively.  4 

We  now  look  at  the  left  boundary  of  a  PZT  pair.  Recall  that  the  beam  model  is 
actually  composed  of  different  PDEs  “pasted”  together  using  boundary  conditions. 
Consider  a  left  boundary  at  x  =  A.  Let  £  be  a  integer  such  that  xi-\  <  A  <  X£. 


4If  the  stiffness  of  the  beam  is  much  greater  than  its  damping,  (i.e.  Eb  E£),  then  the  ODEs 
given  by  Equations  B.80  and  B.83  can  be  approximated  by  the  following  algebraic  expressions 
for  all  (t  >  0): 


uw+i  =  -wN- i  +  2  uiN 
wn+2  =  —  2i0jv-i  +  3wn- 


As  a  result,  the  ODEs  for  the  two  rightmost  points  on  the  beam  become 

(  p  _|^  p*  _d_  \ 

b  bdt  )  {wN- 3  -  4uw-2  +  5wjv_i  -  2u>n) 


Wn-i  =  - 


A4 


wN  =  - 


rh  +  r 


b  dt 


A4 


Ojv-2  -  2u>at-i  +  WN) . 


(B.86) 

(B.87) 

(B.88) 

(B.89) 
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We  assume  that  there  is  no  other  boundary  “nearby”  so  that,  for  points  of  interest 
to  the  boundary,  x%  is  in  a  simple  section  for  i  <  i  and  and  ay  is  in  a  sandwiched 
section  for  i  >  £.  The  the  displacement  expression  given  in  Equation  B.71  is  not 
valid  at  The  simple  section  points  xe-2  and  xe-i  are  boundary  points  since  the 
displacement  expression  given  in  Equation  B.71  requires  the  displacement  values 
at  points  outside  the  section.  Likewise,  the  displacement  expression  given  in  Equa¬ 
tion  B.72  is  not  valid  for  the  sandwiched  section  points  xi  and  xe+\-  Hence,  we 
introduce  ghost  points  X£,  and  xe+i  (to  augment  the  simple  section)  and  xe~2  and 
xt.x  (to  augment  the  sandwiched  section). 

Now  we  use  the  boundary  conditions  given  in  Equations  B.51  through  B.58  to 
derive  expressions  for  the  ghost  displacements  we-2,  wi- 1,  we,  and  We+i  in  terms 
of  the  actual  displacements.  First,  we  introduce  some  notation  to  simplify  the 
resulting  expressions.  Let 


Ebh  +  ECIC 
Ebh 


(B.90) 


and 


A  =  - bd3iEc(tb  +  tc). 


(B.91) 


The  boundary  condition 


w( A  ,t)  =  w(X+,t)  Vt  >  0 


(B.92) 


implies 


tv?  =  we  Vf  >  0. 


The  boundary  condition 


dw 

dx 


(A  ,t) 


dw 

dx 


(AV) 


Vt  >  0 


implies 


we+i  -we_  we+i  -  we 

A  “  A 


(B.93) 


(B.94) 
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=»  w(+i  =  we+ 1  Vf  >  0. 


(B.95) 


The  boundary  condition  given  by  Equation  B.53  becomes 


f)2  yj  f)2  yj 

—  (A-,f)  =  -Af/c  +  T  — (A+,t)  \/t  >  0, 


(B.96) 


which  implies 


we- 1  -  2u>£  +  we+i  =  -A2AVC  +  4/  (r^_i  -  2we  +  icm) 


—  2w(  +  +  A2AVC  +  24/w;^  —  =  4/'i%_1 

A2 A  1  2(4/  —  1)  T  - 1 

=  vr14 + +  —* — ~  -y~Wt+i 


(B.97) 


for  all  t  >  0.  The  boundary  condition  given  by  Equation  B.54  becomes 


!>-•*> =»§>+-*>  vi>«- 


(B.98) 


which  implies 


-w£-2  +  3we-i  -  3we  +  we+ 1  =  4/  (-we-2  +  3w(-i  -  3 we  +  wg+i) 


-we- 2  +  3we-i  -  3u^  +  t^+i  -  3'S'we-i  +  3'S'we  -  d'uy+i  = 


=*  we~2  =  yWe- 2  - 

/  A2A  1  2(4/  —  1)  T  —  1  A 

+3  I  — bc  +  — H - - - we - — we+i  1  -3 we  +  we+ 1 


3A2At,  1  3(4/  —  1)  2(4/  —  1)  , 

=»  We~2  =  —y—Vc  +  ^wt~2  4 - ^ - wt - ^ - we+i-  (B.99) 

Hence,  to  the  left  of  the  boundary  we  get 


*f-  e*a 

we- 2  = - bdt  (we-4  -  4wd-3  +  6u^_2  -  4 wt-e  +  we) 


rb+ni 


(we- 4  -  4tc^_3  +  6w(-2  -  4 we-i  +  we)  (B.100) 
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and 


Wl-2  = 


A4 

r6  +  vti 

A4 


*  d 

bdt  (we- 3  -  4:W£—2  +  6we-i  -  Awe  +  we+i) 
bdt  (we- 3  -  4u^_2  +  6^_i  -  4-u;^  +  t^+i) .  (B.101) 


To  the  right  of  the  boundary 


wt  = 


p  p*if 

c  dt  (we-2  -  Awt-i  +  6iu*  -  4uy+1  +  we+2 ) 


A4 

rc  +  r 


*  d_ 
c  dt 


-4 


A4 

/  A2A 


/  3A2A 


1 


vp  rc  +  ?«.t_2  + 


3(T  -  1)  2(T  -  1) 


-we 


-turn 


1  2(*  - 1)  -  1 
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and 


we+i  = 


v  ±r*4 

1  C  T  i  ,  j , 


cdt 


A4 

rc  +  r*|  / a2 a 


(we-i  -  4w£  +  6wy+i  -  4uy+2  +  u^+3) 


A4 
T  -  1 


T,  1  2(T  —  1) 

V  T  K  +  ^-1+  T  ^ 


:u>m  -  Awe  +  6^+i  -  Awe+2  +  wm 


rc  +  r*i  /i 


C  dt 


A4 


+ 


2(T  -  1) 


-4  U< 


6  -  —  -  ^  we+ 1  -  4^+2  +  (B.103) 


Here  we  have  assumed  that  T* ^Vc  is  small  compared  to  TCVC  and  can  be  neglected. 
In  practice,  Tc  T*  and  Vc  is  smooth,  so  this  assumption  is  reasonable. 

A  similar  procedure  is  followed  to  obtain  ODEs  for  the  boundary  points  at  the 
right  edge  of  a  crystal  pair.  The  resulting  equations  are  stated  in  the  following 
section. 


165 


B.4.1  Finite  Difference  Model  Statement 


Here  we  summarize  the  previous  section  by  stating  the  complete  set  of  ODEs  which 
comprise  the  finite  difference  approximation  of  the  beam  model. 

Define  £{i)  to  be  an  integer  such  that  £r(*)-i  <  @(i)  <  %Ui)-  Likewise,  dehne 
r(i)  to  be  an  integer  such  that  xr^  <  r(i)  <  xr^+i.  additionally,  let  r( 0)  =  2  and 
£(k  +  1  )  =  N. 
simple  sections: 

For  i  such  that  r(j)  <  i  <  £{j  +  1),  j  e  {0,1,2,...,  k}, 

P  p*  d 

Wi  = - b  A4  bdt  (wi- 2  -  1  +  Qwt  -  Awi+i  +  wi+2)  •  (B.104) 

sandwiched  sections: 

For  i  such  that  £{j)  +  2  <  i  <  r(j )  —  2,  j  e  {1,  2,  3, ... ,  k}, 

P  p*  d 

Wi  = - C  ^ACdt  (wi- 2  -  4tCi_i  +  Qwt  -  4wi+i  +  wi+2)  •  (B.105) 

boundary  conditions  (clamped  end): 

p  I  p*_d_ 

= - b  ^4bdt  (Gw i  -  4tc2  +  w3) ,  (B.106) 

p  _l  p*  A. 

w2  = - b  ^4bdt  (~4wi  +  6w2  -  4w3  +  wA) .  (B.107) 

boundary  conditions  (free  end): 

Tb  +  r*A 

WN- 1  = - A4  hdt  (wjv-3  -  4tCjv_2  +  6tCjv_i  -  4wjv  +  wN+1)  (B.108) 

Tb  +  r*A 

wN  = - bdt  (wN_2  ~  4tcjv„i  +  6wn  -  4wN+1  +  WN+2) ,  (B.109) 


where  Wn+ 1  and  Wn+2  satisfy 


d  „ 

JtWN+l 

d  _ 

d^+2 


Ebh 

Wi> 

Eblb 

Wb 


(%- 1  —  2u)at  +  icat+i)  —  rhjv— 1  +  2-ih/v 
(2"ICjv— 1  -  3wat  +  fhjv+2)  “  2tCAT-l  +  3tCAr. 


(B.110) 

(B.lll) 
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boundary  conditions  (left  edge  of  PZT  pair): 

For  i  G  {1, 2,3,...,  k}, 


Wi{i) 


rc  +  r:!  (i  4  /  5(^-1)' 

— U^)-2  -  +  ( 6  +  — * —  )  ,rF0 


2(4/  -  1) 
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TCA 

A24> 


v: 


(B.112) 


rc  +  r:-|  /i  / 2(4/  —  i) 

- A?-  ~  4  )  ro<«) 


T  -  1 
4/ 


'«//(*)+!  —  4lC^(j)+2  +  W£(i)+ 3 


4Ar‘.  (B.113) 


boundary  conditions  (right  edge  of  PZT  pair): 

For  i  G  {1, 2,3,...,  k}, 


kVr(i)  —  1 


F  _L  V*A 
1  c  Lcdt 


A4 
2(4/  —  1) 
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A24> 


^r(i)— 1 

Vc  (B.114) 
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r(i) 


p  p*  d 

1  c  ^  1  cdt 


A4 

5(4/  -  1) 
4/ 


wr(i)~  2  + 


2(T  -  1) 

4> 


4  ^V(i)  — 1 


4  1  \  r.A 

-  ^«V(i)+i  +  ^«>r(i)+2  J  +  (B-H5) 


Given  initial  conditions  Wj(0)  and  Wj(0)  for  i  =  {1,2,...,  A},  w;Ar+i(0),  and 
w;v+2(0),  and  inputs  Vj(t),  j  =  { 1,2,...,  k},  \/t  >  0,  these  equations  give  a  com¬ 
plete  description  of  the  deformation  of  the  beam  for  all  time  t  >  0.  Hence,  we 
have  approximated  the  PDEs  listed  in  Section  B.3.1  with  a  set  of  N  second  order 
ODEs  along  with  two  first  order  ODEs. 


B.5  Simulation  Results 

In  this  section  we  present  results  of  the  simulation  of  a  flexible  beam  based  on 
the  above  model.  The  N  second  order  ODEs  listed  in  Section  B.4.1  were  first 
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transformed  into  state-space  form  to  creating  2 N  first  order  ODEs.  These  were 
then  numerically  integrated  using  the  fourth  order  Runge-Kutta  method. 

For  the  simulation,  we  chose  to  model  an  aluminum  beam  with  a  length  of  1 
meter,  a  width  of  0.05  meters,  and  a  thickness  of  2.5  millimeters.  Four  pairs  of 
PZT  crystals  are  evenly  spaced  along  the  length  of  the  beam. 

Figures  B.8  and  B.9  depict  the  response  of  the  beam  to  a  100  volt  step  applied 
to  the  leftmost  actuator  pair  at  time  t  —  0.  Figure  B.8  shows  the  configuration 
of  the  entire  beam  at  ten  successive  “snapshots”  in  time.  Figure  B.9  plots  the 
position  of  the  tip  of  the  beam  as  a  function  of  time.  Both  figures  compare  well 
with  experimentally  observed  results  for  a  similar  beam  [26]. 
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Figure  B.8:  Snapshots  in  time  of  si 
step  input  to  the  leftmost  actuator 
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